In [2]:
import cv2
from matplotlib import pyplot as plt
import numpy as np
import os
import scipy.signal as sgn
In [ ]:
 
In [3]:
# Ivo's implemetation of Ellipse class

class Ellipse:

    C : np.array #homogeneous conic matrix x'*C*x=0

    def __init__( self, C : np.array ):

        self.C = C
        

    @classmethod
    def estimate_ellipsoid( cls, points : np.array ):
        """
        C = estimate_ellipsoid_2D ( points )
        ----------------------------------

        estimates a homogeneous ellipsoid C
        
        x'*C*x = 0
        
        from a number of point samples 
        (in a least squares sense)
        
        format of points

                [ x^1_1 . . . . .  x^N_1 ] 
                [  .                 .   ]
        points = [  .                 .   ]
                [  .                 .   ]
                [ x^1_n . . . . .  x^N_n ]
                [  1    . . . . .    1   ]

        where each point is element R^n
        and there are N points, the points  
        have to be in homogeneous coordinates.
        
        """

        n = points.shape[0]
        print(n)
        num_entries = int((n+1) * n/2)
        
        if points.shape[1] < num_entries-1:
            print('insufficient number of points, need at least {}'.format( num_entries-1 ) )
    
        # --------------  estimate ellipse ---------------

        # assemble linear system
        A = np.zeros( [points.shape[ 1 ], num_entries ] )

        #for all points
        for i in range( 0, points.shape[ 1 ] ):
        
            #extract point
            x = np.array([points[:, i]]).T
            #print(x)
            X  = np.matmul( x, x.T )
            #print(X)
            X2 = X.T + X - np.diag( np.diag( X ) )

            #print(X2)
            row = X2[ np.where( np.triu( np.ones( X2.shape ) ) ) ]
            #print(row)

            A[ i, : ] = row

        #print(A)

        #solve for null space 
        U,D,V = np.linalg.svd( A )

        ell = V[-1,:]


        #create conic matrix
        Cd = np.zeros( X2.shape )

        Cd[ np.where( np.triu( np.ones( X2.shape ) ) ) ] = ell

        C = Cd + Cd.T - np.diag( np.diag ( Cd ) )

        #detC =  np.linalg.det( C[0:-1,0:-1] )

        #print(detC)

        if Ellipse.is_hyperboloid( C ):
            print('Warning: result is a hyperboloid!')
        

        if Ellipse.is_paraboloid( C ):
            print('Warning: result is a paraboloid!')
        
        return C

    @classmethod
    def cov_mean_to_ellipsoid( cls, cov, mean, k ):
        """
            C = cov_mean_to_ellipsoid_homog ( Cov, Mean, k )
            ----------------------------------
        
            generates a homogeneous ellipsoid C
            
            x'*C*x = 0
        
            from covariance matrix and mean value
            of a gaussian distribution
            
            k is a confidence factor determibed by the chi^2 
            distribution. The probability that x is inside 
            the ellipsoid is P_{chi^2}(k,2). Some values 
            for k are 2.41 (70%) and k=5.99 (95%).
            
            cov  : np.array - covariance matrix (non-homogeneous)
            mean : np.array - mean value (vector)
            k : float - confidence factor, see above 

        """

        #column vector, assume we get a 1D array
        mean = np.array( [mean] ).T

        SIGMAinv = np.linalg.inv( cov )

        lastentry = np.matmul( mean.T, np.matmul( SIGMAinv, mean ) ) - k 
        C = np.vstack( [ np.hstack( [ SIGMAinv, -np.matmul(SIGMAinv, mean)] ), 
                         np.hstack( [ -np.matmul(mean.T,SIGMAinv), lastentry] ) ] )
        
        return C

    @classmethod
    def srt_to_ellipsoid( cls, S, R, T ): 
        """
         C = srt_to_ellipsoid_homog ( S, R, T )
         ----------------------------------
        
         generates a homogeneous ellipsoid C
         
         x'*C*x = 0
        
         from translation, rotation and scale parts as
        
         x'* (T^-1 * R^-1 * S^-1) ( S * R * T ) x = 0
        
         see: Ivo Ihrke: "Some Notes on Ellipses", TechReport, 2004
        """
        #function C = srt_to_ellipsoid_homog ( S, R, T )


        invT = np.linalg.inv( T )
        invR = R.T
        invS = np.linalg.inv( S )

        left  = np.matmul( invT.T, np.matmul( invR.T, invS.T ) )
        right = np.matmul( invS, np.matmul( invR, invT ) )

        C = np.matmul( left, right )
    
        return C



    @classmethod
    def decompose_ellipsoid( cls, C ):
        """
        
         [S,R,T] = decompose_ellipsoid_homog ( C )
         ----------------------------------
        
         decomposes a homogeneous ellipsoid C
         
         x'*C*x = 0
        
         into translation, rotation and scale parts
        
         x'* (T^-1 * R^-1 * S^-1) ( S * R * T ) x = 0
        
         Note that C \neq (T^-1 * R^-1 * S^-1) ( S * R * T )
         in some cases because the orientation of 
         the ellipsoid's major axis is not fixed.
         The function generates a right handed system.

        """

        #print(C)
        C_ih = C[0:-1,0:-1] 
        #print(C_ih)

        t = -np.matmul( np.linalg.pinv( C_ih ), C[0:-1,-1] )

        T = np.eye( C.shape[ 0 ] )
        T[0:-1, -1 ] = t

        U,D,V = np.linalg.svd( C_ih )

        R = np.eye( C.shape[ 0 ] )
        R[0:-1,0:-1] = V


        #S = R' * T' * C * T * R;
        S = np.matmul( R.T, np.matmul( T.T, np.matmul( C, np.matmul( T, R ) ) ) )

        S = S / S[-1,-1]

        S = np.sqrt( np.diag(1/abs(np.diag(S))))

        return S, R, T

    @classmethod
    def is_ellipsoid( cls, C ):
        """
         bool = is_ellipsoid_homog ( C )
         ----------------------------------
        
         determines whether homogeneous conic C
         
         x'*C*x = 0
         
         is an ellipsoid
        """

        #function val = conic_is_ellipsoid_homog ( C )

        detC =  np.linalg.det( C[0:-1,0:-1] )

        return detC > 0


    @classmethod
    def is_hyperboloid( cls, C ):
        """
         bool = is_hyperboloid_homog ( C )
         ----------------------------------
        
         determines whether homogeneous conic C
         
         x'*C*x = 0
         
         is a hyperboloid
        """

        
        detC =  np.linalg.det( C[0:-1,0:-1] )

        return detC < 0

    @classmethod
    def is_paraboloid( cls, C ):
        """
         bool = is_paraboloid_homog ( C )
         ----------------------------------
        
         determines whether homogeneous conic C
         
         x'*C*x = 0
         
         is a paraboloid
        """

        
        detC =  np.linalg.det( C[0:-1,0:-1] )

        #maybe with EPS ?
        return detC == 0


    def get_center( self ):
        
        S,R,T = self.decompose_ellipsoid( self.C )


        return [ T[2,0], T[2,1] ]


    def plot( self, ax, **kwargs ):
        """
        plot_ellipsoid_2D ( C, style, N )
        ----------------------------------
        
        plots a homogeneous ellipsoid C
        
        x'*C*x = 0
        
        style is a matlab plot style, e.g. 'b-' (default)
        N is the number of polygon segments to 
        approximate the ellipse
        """

        #function plot_ellipsoid_homog_2d ( C, style, N )

        if not 'style' in locals():
            style = 'b-'

        if not 'N' in locals():
            N = 100



        angles = np.linspace( -np.pi, np.pi, N )

        # generate a unit circle or a unit hyperbola
        if ( Ellipse.is_ellipsoid( self.C ) ):
            #uc  = [ cos(x); sin(x); ones(size(x)) ];
            x  = np.cos( angles )
            y  = np.sin( angles )
            w  = np.ones( angles.shape )

            uc = np.vstack([x,y,w])

        else:
            EPS = 0.1
            x   =  np.linspace( -np.pi/2+EPS, np.pi/2-EPS, int( N/2 ) )#-(pi/2-EPS):((pi-2*EPS)/(N/2)):(pi/2-EPS);
            x2  =  np.linspace( np.pi/2+EPS, 3*np.pi/2-EPS, int( N/2 ) )  #x2 =  (pi/2+EPS):((pi-2*EPS)/(N/2)):(3*pi/2-EPS);
            uc  = np.vstack( [ 1/np.cos(x), np.tan(x), np.ones( x.shape ) ] )
            uc2 = np.vstack( [ 1/np.cos(x2), np.tan(x2), np.ones( x2.shape ) ] )
        
        #decompose
        S,R,T = Ellipse.decompose_ellipsoid( self.C )

        #transform unit circle or parabola 
        transform = np.matmul( T, np.matmul( R, S ) )
        ell  =  np.matmul( transform, uc )

        if 'uc2' in locals():
            ell2  =  np.matmul( transform, uc2 )
            ax.plot( ell2[0,:], ell2[1,:], style )

   
        #plot conic
        ax.plot( ell[0,:], ell[1,:], style )

        #plot coordinate system
        tx = T[0,2]
        ty = T[1,2]
        
#         print(tx, ty)

        smin = S[0,0]
        smaj = S[1,1] 

        ax.plot( np.array( [ tx, tx+smin*R[0,0] ] ), np.array( [ ty, ty+smin*R[1,0] ] ), 'r' )
        ax.plot( np.array( [ tx, tx+smaj*R[0,1] ] ), np.array( [ ty, ty+smaj*R[1,1] ] ), 'b' )

        ax.set_aspect('equal', adjustable='box')
        return (tx, ty)
In [4]:
def fit_ellipse(path_dir_w_time, ifHull=False, date_now=""):
        """ Find and fit an ellipse in an image.
        it should be run on a normalized image I (normalization should exclude dead pixels which are artificially maximum)
        Parameters
        ----------
        path_dir_w_time : str
            _description_
        date_now : str
            _description_

        Returns
        -------
        Ellipse
            Ellipse in an image.
        """
        
        I_temp = plt.imread(path_dir_w_time)[:,:,1].astype(float)/255
        I = old_rplc_dead_px_med(I_temp)

        ellipsoid_points=[]

        if ifHull == False:
            #every Nth scanline should be sampled
            START_SCANLINE = 50 #there are some weird effects on the zero scanline that lead to detection of erroneous points
            DELTA_SCANLINE = 100

            #need to filter enough such that the big jump into and out of the disc has indeed the largest magnitude derivative
            FILTERWIDTH=250
            x=np.linspace(-4,4,FILTERWIDTH)
            gauss=np.exp(-x**2)


            plt.figure(figsize=[18,12])

            for sl in range( START_SCANLINE, I.shape[0]-START_SCANLINE, DELTA_SCANLINE ):
                #print(f'Checking Scanline {sl} ...' )
                #line = self.I[sl]
                line = I[sl, :]

                filtered_line = sgn.convolve(line,gauss,mode='same')
                #crude filter for line intersecting the ring structure
                if np.amax(filtered_line) > 0.9:
                    #print("- scanline intersection ")

                    filtered_deriv = np.diff( filtered_line )

                    #crude extraction - this will probably need improvement
                    max_pos = np.argmax( filtered_deriv )
                    min_pos = np.argmin( filtered_deriv )
                    ellipsoid_points.append([max_pos,sl])
                    ellipsoid_points.append([min_pos,sl]) #add two points on perimeter

                else:
                    #print("- no intersection ")

                    pass
            print(ellipsoid_points)


        else:
            img = np.array(I*255, dtype=np.uint8)
            cnt_img, final_cnt = contour_points(img)

            if final_cnt is None:
                return None, None
            elif final_cnt.shape[0] < 5:
                return None, None
            else:
                ellipsoid_points = np.array(final_cnt).reshape([-1, 2])
        
#         print(ellipsoid_points)
        ellipsoid_points = np.array(ellipsoid_points).astype('float').T
        ellipsoid_points = np.vstack( [ellipsoid_points, np.ones(ellipsoid_points.shape[1])] )
        print(ellipsoid_points)

        #plt.clf()
        #plt.imshow(self.I)
        #plt.show()

        for i in range( 0, ellipsoid_points.shape[1] ):
            plt.plot(ellipsoid_points[0,i],ellipsoid_points[1,i],'rx')

        E = Ellipse( Ellipse.estimate_ellipsoid( ellipsoid_points ) )
        S,R,T=Ellipse.decompose_ellipsoid(E.C)

        plt.plot(T[0,2],T[1,2],'gx')

        plt.imshow(I)
        center = E.plot(plt.gca())
        plt.title(f'ellipsoide fit {1}')
#         plt.savefig(f"{path_dir_w_time}/{date_now}_ledpos_{self.x_pos}_{self.y_pos}_{self.channel}_ellipsoide_fit.png")

        return E, center

def old_rplc_dead_px_med(I_temp):
        """Replace dead pixels by median

        Parameters
        ----------
        I_temp : ndarray
            image to be processed

        Returns
        -------
        ndarray
            clean image where dead pixels are replaced by median 
        """

        Imed = sgn.medfilt2d( I_temp, kernel_size=5 )
        maxI = np.amax( Imed ) #normalize robustly
        I_temp = I_temp / maxI
        #replace dead pixels by median filtered values
        I_temp[ I_temp > 1 ] = Imed[ I_temp > 1 ]

        return I_temp
In [5]:
def contour_points(g_img):

    kernel = np.ones((3, 3), np.uint8)
    gradient = cv2.morphologyEx(g_img, cv2.MORPH_GRADIENT, kernel)

    blur = cv2.GaussianBlur(gradient, (3, 3), 0)
    ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)

    gradient_img = np.zeros_like(g_img, dtype=np.uint8)

    contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
    if len(contours) < 2:
        return None, None
    
    cv2.drawContours(gradient_img, contours, -1, 255, 2)

    kernel = np.ones((5, 5), np.uint8)
    cv2.dilate(gradient_img, kernel, iterations = 1)


    hull_img = np.zeros_like(g_img, dtype=np.uint8)

    blur = cv2.GaussianBlur(g_img, (3, 3), 0)
    ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)

    contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
    
    if len(contours) < 2:
        return None, None
    
    area_cntrs = [cv2.contourArea(cnt) for cnt in contours]
    max_cntr_index = np.argmax(area_cntrs)
    cnt = contours[max_cntr_index]

    hull = cv2.convexHull(cnt)
    cv2.drawContours(hull_img, [hull], 0, 125, 2)

    op = np.logical_and(gradient_img, hull_img).astype(np.uint8)*255

    op[:10, :] = 0
    op[-10:, :] = 0

    kernel = np.ones((3, 3), np.uint8)
    opening = cv2.morphologyEx(op, cv2.MORPH_OPEN, kernel)
    
    x, y = np.where(opening > 1)
    final_cnt = np.vstack([y, x]).T
#     final_cnt = np.take(final_cnt, list(set(np.random.randint(0, final_cnt.shape[0], size=int(final_cnt.shape[0]/10)))), axis=0)
    imgs = np.hstack([gradient_img, hull_img, opening])
    
    if final_cnt.shape[0] > 4:
        return opening, final_cnt
    
    else:
        return None, None
In [ ]:
 

Problem with Preprocessing: Demo¶

In [233]:
dir_path = "Z:/CSE/CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_08_25/2023_08_25_09_31_25/"
file_path = "2023_08_25_09_31_25_img_iso_100_shutter_2_x_16_y_16_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img, _ = cv2.split(img)
In [238]:
imgs, cnt = contour_points(g_img)

kernel = np.ones((25, 25),np.uint8)
cv2.dilate(imgs, kernel, iterations = 1)

fig, ax = plt.subplots(1, 1, figsize=(15, 15))
plt.imshow(imgs, cmap='gray')
plt.show()

ellipse = cv2.fitEllipse(cnt)
cv2.ellipse(g_img, ellipse, 255, 2)
fig, ax = plt.subplots(1, 1, figsize=(15, 15))
ax.imshow(g_img, cmap='gray')
plt.show()

print(ellipse)
((2971.643798828125, 913.9298095703125), (5347.255859375, 5739.8564453125), 169.41822814941406)
In [ ]:
 
In [ ]:
 

Preprocessing problem fixed: Demo¶

In [652]:
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_13_y_14_r_0_g_1_b_0.tiff"

# dir_path = "Z:/CSE/CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_08_25/2023_08_25_09_31_25/"
# file_path = "2023_08_25_09_31_25_img_iso_100_shutter_2_x_16_y_16_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)

_, g_img, _ = cv2.split(img)
In [653]:
g_img
Out[653]:
array([[   0,    0,    0, ...,    0,    0,    0],
       [   0,    0,    0, ...,    0,    0,    0],
       [   0,    0,    0, ...,    0,    0,    0],
       ...,
       [ 887,  153, 1258, ...,    0,    0,    0],
       [1531, 1343, 2342, ...,    0,    0,    0],
       [1369, 1267, 1454, ...,    0,    0,    0]], dtype=uint16)
In [654]:
plt.imshow(g_img, cmap='gray')
Out[654]:
<matplotlib.image.AxesImage at 0x208ad516f10>
In [655]:
g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
In [656]:
plt.imshow(g_img_8, cmap='gray')
Out[656]:
<matplotlib.image.AxesImage at 0x208b074b070>
In [657]:
kernel = np.ones((15, 15), np.uint8)
dilation = cv2.dilate(g_img_8, kernel, iterations = 1)
plt.imshow(dilation, cmap='gray')
plt.show()

kernel = np.ones((15, 15), np.uint8)
erosion = cv2.erode(g_img_8, kernel, iterations = 1)
plt.imshow(erosion, cmap='gray')
plt.show()

plt.imshow((np.fabs(-erosion.astype(np.int16) + dilation.astype(np.int16))).astype(np.uint8), cmap='gray')
plt.show()
In [ ]:
 
In [658]:
kernel = np.ones((15, 15), np.uint8)
gradient = cv2.morphologyEx(g_img_8, cv2.MORPH_GRADIENT, kernel)

# g_img_2 = g_img.cop
plt.imshow(gradient, cmap='gray')
Out[658]:
<matplotlib.image.AxesImage at 0x208b23ed430>
In [659]:
# kernel = np.ones((25, 25), np.uint8)
# opening = cv2.morphologyEx(gradient, cv2.MORPH_OPEN, kernel)

# plt.imshow(opening, cmap='gray')
In [668]:
blur = cv2.GaussianBlur(gradient+20, (5, 5), 0)
ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)

gradient_img = np.zeros_like(g_img, dtype=np.uint8)

contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
if len(contours) < 2:
    print("Errors")
    
cv2.drawContours(gradient_img, contours, -1, 255, 2)

kernel = np.ones((35, 35),np.uint8)
gradient_img = cv2.dilate(gradient_img, kernel, iterations = 1)

plt.imshow(gradient_img, cmap='gray')
Out[668]:
<matplotlib.image.AxesImage at 0x208d1f2adf0>
In [669]:
ch_img_1 = np.zeros_like(img, dtype=np.uint8)
ch_img_1[:, :, 0] = gradient_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_8+50

new_img = np.array((ch_img_1 + ch_img_2), dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [ ]:
 
In [670]:
hull_img = np.zeros_like(g_img_8, dtype=np.uint8)

blur = cv2.GaussianBlur(g_img_8, (3, 3), 0)
ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)

plt.imshow(th3, cmap='gray')
Out[670]:
<matplotlib.image.AxesImage at 0x208d0ee2070>
In [ ]:
 
In [671]:
contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
area_cntrs = [cv2.contourArea(cnt) for cnt in contours]
max_cntr_index = np.argmax(area_cntrs)
cnt = contours[max_cntr_index]
hull = cv2.convexHull(cnt)


if len(contours) < 2:
    print("Errors")
    
cv2.drawContours(hull_img, [hull], 0, 125, 2)

plt.imshow(hull_img, cmap='gray')
Out[671]:
<matplotlib.image.AxesImage at 0x208d0f92f10>
In [672]:
# hull_img = np.zeros_like(g_img_8, dtype=np.uint8)

# area_cntrs = [cv2.contourArea(cnt) for cnt in contours]
# max_cntr_index = np.argmax(area_cntrs)
# cnt = contours[max_cntr_index]

# hull = cv2.convexHull(cnt)
# cv2.drawContours(hull_img, [hull], 0, 125, 2)
# plt.imshow(hull_img, cmap='gray')
In [673]:
ch_img_1 = np.zeros_like(img, dtype=np.uint8)
ch_img_1[:, :, 0] = hull_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_8

new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [674]:
op = np.logical_and(gradient_img, hull_img).astype(np.uint8)*255

op[:10, :] = 0
op[-10:, :] = 0

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(op, cmap='gray')
plt.show()
In [675]:
# op = np.logical_and(gradient_img, hull_img).astype(np.uint8)*255

# op[:10, :] = 0
# op[-10:, :] = 0

ch_img_1 = np.zeros_like(img, dtype=np.uint8)
ch_img_1[:, :, 0] = op

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_8+50

new_img = np.array((ch_img_1 + ch_img_2), dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [228]:
# ellipse fit below not terrific but not because of contour above is v. good
In [676]:
x, y = np.where(op > 1)
op_cnt = np.vstack([y, x]).T

ellipse = cv2.fitEllipse(op_cnt)
cv2.ellipse(img, ellipse, (2**15, 0, 0), 3)


fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(cv2.convertScaleAbs(img, alpha=(255/65535)) + 125)
plt.show()
In [230]:
ellipse
Out[230]:
((2975.01220703125, 1009.5517578125),
 (5253.75927734375, 5340.55322265625),
 97.18701934814453)
In [231]:
min_axis_len = 5253.75927734375
maj_axis_len = 5340.55322265625

ratio = min_axis_len/maj_axis_len
ratio
Out[231]:
0.9837481358777039
In [ ]:
 
In [6]:
def contour_points_v2(g_img):

    kernel = np.ones((15, 15), np.uint8)
    gradient = cv2.morphologyEx(g_img, cv2.MORPH_GRADIENT, kernel)
    
    blur = cv2.GaussianBlur(gradient + 20, (5, 5), 0)
    ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)

    gradient_img = np.zeros_like(g_img, dtype=np.uint8)

    contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
    if len(contours) < 1:
        return None, None
    
    cv2.drawContours(gradient_img, contours, -1, 255, 2)

    kernel = np.ones((25, 25),np.uint8)
    gradient_img = cv2.dilate(gradient_img, kernel, iterations = 1)    


    hull_img = np.zeros_like(g_img, dtype=np.uint8)

    blur = cv2.GaussianBlur(g_img, (3, 3), 0)
    ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)

    contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
    
    if len(contours) < 1:
        return None, None
    
    area_cntrs = [cv2.contourArea(cnt) for cnt in contours]
    max_cntr_index = np.argmax(area_cntrs)
    cnt = contours[max_cntr_index]

    hull = cv2.convexHull(cnt)
    cv2.drawContours(hull_img, [hull], 0, 125, 2)

    op = np.logical_and(gradient_img, hull_img).astype(np.uint8)*255

#     op[:10, :] = 0
#     op[-10:, :] = 0

#     kernel = np.ones((3, 3), np.uint8)
#     opening = cv2.morphologyEx(op, cv2.MORPH_OPEN, kernel)
    
    x, y = np.where(op > 1)
    final_cnt = np.vstack([y, x]).T
#     final_cnt = np.take(final_cnt, list(set(np.random.randint(0, final_cnt.shape[0], size=int(final_cnt.shape[0]/10)))), axis=0)
#     imgs = np.hstack([gradient_img, hull_img, opening])
    
    if final_cnt.shape[0] > 4:
        return op, final_cnt
    
    else:
        return None, None
In [239]:
# now preprocessing fixed for both types of images
In [ ]:
 
In [ ]:
 

Fixing ellipses' minor and major axis ratio¶

In [532]:
# Ivo's implemetation of Ellipse class

class Ellipse_Fixed_Ratio:

    C : np.array #homogeneous conic matrix x'*C*x=0

    def __init__( self, C : np.array ):

        self.C = C
        

    @classmethod
    def estimate_ellipsoid( cls, points : np.array ):
        """
        C = estimate_ellipsoid_2D ( points )
        ----------------------------------

        estimates a homogeneous ellipsoid C
        
        x'*C*x = 0
        
        from a number of point samples 
        (in a least squares sense)
        
        format of points

                [ x^1_1 . . . . .  x^N_1 ] 
                [  .                 .   ]
        points = [  .                 .   ]
                [  .                 .   ]
                [ x^1_n . . . . .  x^N_n ]
                [  1    . . . . .    1   ]

        where each point is element R^n
        and there are N points, the points  
        have to be in homogeneous coordinates.
        
        """

        n = points.shape[0]
#         n is 3
        num_entries = int((n+1) * n/2) - 2 # fixing ellipse maj/min axis length ratio
        
        if points.shape[1] < num_entries-1:
            print('insufficient number of points, need at least {}'.format( num_entries-1 ) )
    
        # --------------  estimate ellipse ---------------

        # assemble linear system
        A = np.zeros( [points.shape[ 1 ], num_entries ] )

        #for all points
#         print(points.shape[ 1 ])
        for i in range( 0, points.shape[ 1 ] ):
        
            #extract point
            x = np.array([points[:, i]]).T

            X  = np.matmul( x, x.T )

            X2 = X.T + X - np.diag( np.diag( X ) )

            row = X2[ np.where( np.triu( np.ones( X2.shape ) ) ) ]

            

            row[0] = row[0] + row[3]
            row = np.delete(row, 3)
            row = np.delete(row, 1)

            A[ i, : ] = row

        #solve for null space 
        U,D,V = np.linalg.svd( A )
#         print(U, D, V)

        ell = V[-1,:]
        ell = np.insert(ell, 1, 0)
        ell = np.insert(ell, 3, ell[0])
#         ell[0] = 0.9837*ell[0]
        print('Ell:', ell)

        #create conic matrix
        Cd = np.zeros( X2.shape )

        Cd[ np.where( np.triu( np.ones( X2.shape ) ) ) ] = ell

        C = Cd + Cd.T - np.diag( np.diag ( Cd ) )
        print(f'C: {C}')

        if Ellipse_Fixed_Ratio.is_hyperboloid( C ):
            print('Warning: result is a hyperboloid!')
        

        if Ellipse_Fixed_Ratio.is_paraboloid( C ):
            print('Warning: result is a paraboloid!')
        
        return C

    @classmethod
    def cov_mean_to_ellipsoid( cls, cov, mean, k ):
        """
            C = cov_mean_to_ellipsoid_homog ( Cov, Mean, k )
            ----------------------------------
        
            generates a homogeneous ellipsoid C
            
            x'*C*x = 0
        
            from covariance matrix and mean value
            of a gaussian distribution
            
            k is a confidence factor determibed by the chi^2 
            distribution. The probability that x is inside 
            the ellipsoid is P_{chi^2}(k,2). Some values 
            for k are 2.41 (70%) and k=5.99 (95%).
            
            cov  : np.array - covariance matrix (non-homogeneous)
            mean : np.array - mean value (vector)
            k : float - confidence factor, see above 

        """

        #column vector, assume we get a 1D array
        mean = np.array( [mean] ).T

        SIGMAinv = np.linalg.inv( cov )

        lastentry = np.matmul( mean.T, np.matmul( SIGMAinv, mean ) ) - k 
        C = np.vstack( [ np.hstack( [ SIGMAinv, -np.matmul(SIGMAinv, mean)] ), 
                         np.hstack( [ -np.matmul(mean.T,SIGMAinv), lastentry] ) ] )
        
        return C

    @classmethod
    def srt_to_ellipsoid( cls, S, R, T ): 
        """
         C = srt_to_ellipsoid_homog ( S, R, T )
         ----------------------------------
        
         generates a homogeneous ellipsoid C
         
         x'*C*x = 0
        
         from translation, rotation and scale parts as
        
         x'* (T^-1 * R^-1 * S^-1) ( S * R * T ) x = 0
        
         see: Ivo Ihrke: "Some Notes on Ellipses", TechReport, 2004
        """
        #function C = srt_to_ellipsoid_homog ( S, R, T )


        invT = np.linalg.inv( T )
        invR = R.T
        invS = np.linalg.inv( S )

        left  = np.matmul( invT.T, np.matmul( invR.T, invS.T ) )
        right = np.matmul( invS, np.matmul( invR, invT ) )

        C = np.matmul( left, right )
    
        return C



    @classmethod
    def decompose_ellipsoid( cls, C ):
        """
        
         [S,R,T] = decompose_ellipsoid_homog ( C )
         ----------------------------------
        
         decomposes a homogeneous ellipsoid C
         
         x'*C*x = 0
        
         into translation, rotation and scale parts
        
         x'* (T^-1 * R^-1 * S^-1) ( S * R * T ) x = 0
        
         Note that C \neq (T^-1 * R^-1 * S^-1) ( S * R * T )
         in some cases because the orientation of 
         the ellipsoid's major axis is not fixed.
         The function generates a right handed system.

        """

        #print(C)
        C_ih = C[0:-1,0:-1] 
        #print(C_ih)

        t = -np.matmul( np.linalg.pinv( C_ih ), C[0:-1,-1] )

        T = np.eye( C.shape[ 0 ] )
        T[0:-1, -1 ] = t

        U,D,V = np.linalg.svd( C_ih )

        R = np.eye( C.shape[ 0 ] )
        R[0:-1,0:-1] = V


        #S = R' * T' * C * T * R;
        S = np.matmul( R.T, np.matmul( T.T, np.matmul( C, np.matmul( T, R ) ) ) )

        S = S / S[-1,-1]

        S = np.sqrt( np.diag(1/abs(np.diag(S))))

        return S, R, T

    @classmethod
    def is_ellipsoid( cls, C ):
        """
         bool = is_ellipsoid_homog ( C )
         ----------------------------------
        
         determines whether homogeneous conic C
         
         x'*C*x = 0
         
         is an ellipsoid
        """

        #function val = conic_is_ellipsoid_homog ( C )

        detC =  np.linalg.det( C[0:-1,0:-1] )

        return detC > 0


    @classmethod
    def is_hyperboloid( cls, C ):
        """
         bool = is_hyperboloid_homog ( C )
         ----------------------------------
        
         determines whether homogeneous conic C
         
         x'*C*x = 0
         
         is a hyperboloid
        """

        
        detC =  np.linalg.det( C[0:-1,0:-1] )

        return detC < 0

    @classmethod
    def is_paraboloid( cls, C ):
        """
         bool = is_paraboloid_homog ( C )
         ----------------------------------
        
         determines whether homogeneous conic C
         
         x'*C*x = 0
         
         is a paraboloid
        """

        
        detC =  np.linalg.det( C[0:-1,0:-1] )

        #maybe with EPS ?
        return detC == 0


    def get_center( self ):
        
        S,R,T = self.decompose_ellipsoid( self.C )


        return [ T[2,0], T[2,1] ]


    def plot( self, ax, **kwargs ):
        """
        plot_ellipsoid_2D ( C, style, N )
        ----------------------------------
        
        plots a homogeneous ellipsoid C
        
        x'*C*x = 0
        
        style is a matlab plot style, e.g. 'b-' (default)
        N is the number of polygon segments to 
        approximate the ellipse
        """

        #function plot_ellipsoid_homog_2d ( C, style, N )

        if not 'style' in locals():
            style = 'b-'

        if not 'N' in locals():
            N = 100



        angles = np.linspace( -np.pi, np.pi, N )

        # generate a unit circle or a unit hyperbola
        if ( Ellipse_Fixed_Ratio.is_ellipsoid( self.C ) ):
            #uc  = [ cos(x); sin(x); ones(size(x)) ];
            x  = np.cos( angles )
            y  = np.sin( angles )
            w  = np.ones( angles.shape )

            uc = np.vstack([x,y,w])

        else:
            EPS = 0.1
            x   =  np.linspace( -np.pi/2+EPS, np.pi/2-EPS, int( N/2 ) )#-(pi/2-EPS):((pi-2*EPS)/(N/2)):(pi/2-EPS);
            x2  =  np.linspace( np.pi/2+EPS, 3*np.pi/2-EPS, int( N/2 ) )  #x2 =  (pi/2+EPS):((pi-2*EPS)/(N/2)):(3*pi/2-EPS);
            uc  = np.vstack( [ 1/np.cos(x), np.tan(x), np.ones( x.shape ) ] )
            uc2 = np.vstack( [ 1/np.cos(x2), np.tan(x2), np.ones( x2.shape ) ] )
        
        #decompose
        S,R,T = Ellipse_Fixed_Ratio.decompose_ellipsoid( self.C )

        #transform unit circle or parabola 
        transform = np.matmul( T, np.matmul( R, S ) )
        ell  =  np.matmul( transform, uc )

        if 'uc2' in locals():
            ell2  =  np.matmul( transform, uc2 )
            ax.plot( ell2[0,:], ell2[1,:], style )

   
        #plot conic
        ax.plot( ell[0,:], ell[1,:], style )

        #plot coordinate system
        tx = T[0,2]
        ty = T[1,2]
        
#         print(tx, ty)

        smin = S[0,0]
        smaj = S[1,1] 

        ax.plot( np.array( [ tx, tx+smin*R[0,0] ] ), np.array( [ ty, ty+smin*R[1,0] ] ), 'r' )
        ax.plot( np.array( [ tx, tx+smaj*R[0,1] ] ), np.array( [ ty, ty+smaj*R[1,1] ] ), 'b' )

        ax.set_aspect('equal', adjustable='box')
        return (tx, ty, smin, smaj)
In [7]:
def fit_ellipse_2(path_dir_w_time, ifHull=False, date_now=""):
        """ Find and fit an ellipse in an image.
        it should be run on a normalized image I (normalization should exclude dead pixels which are artificially maximum)
        Parameters
        ----------
        path_dir_w_time : str
            _description_
        date_now : str
            _description_

        Returns
        -------
        Ellipse
            Ellipse in an image.
        """
        
        I_temp = plt.imread(path_dir_w_time)[:,:,1].astype(float)/255
        I = old_rplc_dead_px_med(I_temp)

        ellipsoid_points=[]

        if ifHull == False:
            #every Nth scanline should be sampled
            START_SCANLINE = 50 #there are some weird effects on the zero scanline that lead to detection of erroneous points
            DELTA_SCANLINE = 100

            #need to filter enough such that the big jump into and out of the disc has indeed the largest magnitude derivative
            FILTERWIDTH=250
            x=np.linspace(-4,4,FILTERWIDTH)
            gauss=np.exp(-x**2)


            plt.figure(figsize=[18,12])

            for sl in range( START_SCANLINE, I.shape[0]-START_SCANLINE, DELTA_SCANLINE ):
                #print(f'Checking Scanline {sl} ...' )
                #line = self.I[sl]
                line = I[sl, :]

                filtered_line = sgn.convolve(line,gauss,mode='same')
                #crude filter for line intersecting the ring structure
                if np.amax(filtered_line) > 0.9:
                    #print("- scanline intersection ")

                    filtered_deriv = np.diff( filtered_line )

                    #crude extraction - this will probably need improvement
                    max_pos = np.argmax( filtered_deriv )
                    min_pos = np.argmin( filtered_deriv )
                    ellipsoid_points.append([max_pos,sl])
                    ellipsoid_points.append([min_pos,sl]) #add two points on perimeter

                else:
                    #print("- no intersection ")

                    pass
            print(ellipsoid_points)


        else:
            img = np.array(I*255, dtype=np.uint8)
            cnt_img, final_cnt = contour_points_v2(img)

            if final_cnt is None:
                return None, None
            elif final_cnt.shape[0] < 5:
                return None, None
            else:
                ellipsoid_points = np.array(final_cnt).reshape([-1, 2])
        
#         print(ellipsoid_points)
        ellipsoid_points = np.array(ellipsoid_points).astype('float').T
        ellipsoid_points = np.vstack( [ellipsoid_points, np.ones(ellipsoid_points.shape[1])] )
#         print(ellipsoid_points)

        #plt.clf()
        #plt.imshow(self.I)
        #plt.show()

        for i in range( 0, ellipsoid_points.shape[1] ):
            plt.plot(ellipsoid_points[0,i],ellipsoid_points[1,i],'rx')

        E = Ellipse_Fixed_Ratio( Ellipse_Fixed_Ratio.estimate_ellipsoid( ellipsoid_points ) )
        S,R,T=Ellipse_Fixed_Ratio.decompose_ellipsoid(E.C)
        print(f'S:{S}, \nR: {R}, \nT: {T}')

        plt.plot(T[0,2],T[1,2],'gx')

        plt.imshow(I)
        center = E.plot(plt.gca())
        plt.title(f'ellipsoide fit {1}')
#         plt.savefig(f"{path_dir_w_time}/{date_now}_ledpos_{self.x_pos}_{self.y_pos}_{self.channel}_ellipsoide_fit.png")

        return E, center

def old_rplc_dead_px_med(I_temp):
        """Replace dead pixels by median

        Parameters
        ----------
        I_temp : ndarray
            image to be processed

        Returns
        -------
        ndarray
            clean image where dead pixels are replaced by median 
        """

        Imed = sgn.medfilt2d( I_temp, kernel_size=5 )
        maxI = np.amax( Imed ) #normalize robustly
        I_temp = I_temp / maxI
        #replace dead pixels by median filtered values
        I_temp[ I_temp > 1 ] = Imed[ I_temp > 1 ]

        return I_temp
In [257]:
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = f"2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff"

# dir_path = "Z:/CSE/CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_08_25/2023_08_25_09_31_25/"
# file_path = "2023_08_25_09_31_25_img_iso_100_shutter_2_x_15_y_15_r_0_g_1_b_0.tiff"
In [258]:
E, center = fit_ellipse(dir_path + file_path, ifHull=True)
plt.show()

print(center)
[[1.440e+02 1.450e+02 1.460e+02 1.440e+02 1.450e+02 1.460e+02 1.440e+02
  1.450e+02 1.460e+02 4.557e+03 4.558e+03 4.559e+03 4.556e+03 4.557e+03
  4.558e+03 4.559e+03 4.556e+03 4.557e+03 4.558e+03 4.559e+03 4.556e+03
  4.557e+03 4.558e+03 4.480e+03 4.481e+03 4.482e+03 4.480e+03 4.481e+03
  4.482e+03 4.480e+03 4.481e+03 4.482e+03 4.325e+03 4.326e+03 4.327e+03
  4.325e+03 4.326e+03 4.327e+03 4.325e+03 4.326e+03 4.327e+03 3.509e+03
  3.510e+03 3.511e+03 3.501e+03 3.502e+03 3.503e+03 3.504e+03 3.505e+03
  3.509e+03 3.510e+03 3.511e+03 3.496e+03 3.497e+03 3.498e+03 3.499e+03
  3.500e+03 3.501e+03 3.502e+03 3.503e+03 3.504e+03 3.505e+03 3.509e+03
  3.510e+03 3.511e+03 3.492e+03 3.493e+03 3.494e+03 3.495e+03 3.496e+03
  3.497e+03 3.498e+03 3.499e+03 3.500e+03 3.501e+03 3.502e+03 3.503e+03
  3.504e+03 3.505e+03 3.486e+03 3.487e+03 3.488e+03 3.492e+03 3.493e+03
  3.494e+03 3.495e+03 3.496e+03 3.497e+03 3.498e+03 3.499e+03 3.500e+03
  3.486e+03 3.487e+03 3.488e+03 3.492e+03 3.493e+03 3.494e+03 3.495e+03
  3.486e+03 3.487e+03 3.488e+03]
 [1.791e+03 1.791e+03 1.791e+03 1.792e+03 1.792e+03 1.792e+03 1.793e+03
  1.793e+03 1.793e+03 3.192e+03 3.192e+03 3.192e+03 3.193e+03 3.193e+03
  3.193e+03 3.193e+03 3.194e+03 3.194e+03 3.194e+03 3.194e+03 3.195e+03
  3.195e+03 3.195e+03 3.265e+03 3.265e+03 3.265e+03 3.266e+03 3.266e+03
  3.266e+03 3.267e+03 3.267e+03 3.267e+03 3.367e+03 3.367e+03 3.367e+03
  3.368e+03 3.368e+03 3.368e+03 3.369e+03 3.369e+03 3.369e+03 3.730e+03
  3.730e+03 3.730e+03 3.731e+03 3.731e+03 3.731e+03 3.731e+03 3.731e+03
  3.731e+03 3.731e+03 3.731e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03
  3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03 3.732e+03
  3.732e+03 3.732e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03
  3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03 3.733e+03
  3.733e+03 3.733e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03
  3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03 3.734e+03
  3.735e+03 3.735e+03 3.735e+03 3.735e+03 3.735e+03 3.735e+03 3.735e+03
  3.736e+03 3.736e+03 3.736e+03]
 [1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
  1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
  1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
  1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
  1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
  1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
  1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
  1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
  1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
  1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
  1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
  1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
  1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
  1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
  1.000e+00 1.000e+00 1.000e+00]]
3
Warning: result is a hyperboloid!
(5272.143308579376, 3112.6301639181233)
In [259]:
E, center = fit_ellipse_2(dir_path + file_path, ifHull=True)
plt.show()

print(center)
Ell: [ 5.23202987e-07  0.00000000e+00 -1.44310664e-03  5.14674778e-07
 -6.28543710e-04  9.99998761e-01]
C: [[ 5.23202987e-07  0.00000000e+00 -1.44310664e-03]
 [ 0.00000000e+00  5.14674778e-07 -6.28543710e-04]
 [-1.44310664e-03 -6.28543710e-04  9.99998761e-01]]
S:[[2.67648633e+03 0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 2.69857004e+03 0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 1.00000000e+00]], 
R: [[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]], 
T: [[1.00000000e+00 0.00000000e+00 2.75821559e+03]
 [0.00000000e+00 1.00000000e+00 1.22124444e+03]
 [0.00000000e+00 0.00000000e+00 1.00000000e+00]]
(2758.2155895125884, 1221.2444375066455, 2676.486329255531, 2698.5700358857857)
In [260]:
center[2]/center[3]
Out[260]:
0.9918165152889924
In [ ]:
 

Compare multiple images to confirm if there is geometric distrotion¶

In [331]:
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_13_y_15_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img, _ = cv2.split(img)

dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_12_y_15_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img_2, _ = cv2.split(img)

ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [336]:
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + np.roll(ch_img_2, 195, axis=1))*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(20, 20))
ax.imshow(new_img)
plt.show()
In [ ]:
 
In [337]:
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img, _ = cv2.split(img)

dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_14_y_15_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img_2, _ = cv2.split(img)

ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [340]:
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + np.roll(ch_img_2, (195, 0), axis=(1, 0)))*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [ ]:
 
In [346]:
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_17_y_15_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img, _ = cv2.split(img)

dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_16_y_15_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img_2, _ = cv2.split(img)

ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [347]:
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + np.roll(ch_img_2, (195, 0), axis=(1, 0)))*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [ ]:
 
In [348]:
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_19_y_15_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img, _ = cv2.split(img)

dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_18_y_15_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img_2, _ = cv2.split(img)

ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [349]:
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + np.roll(ch_img_2, (195, 0), axis=(1, 0)))*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [ ]:
 
In [350]:
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_13_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img, _ = cv2.split(img)

dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_12_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img_2, _ = cv2.split(img)

ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [355]:
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + np.roll(ch_img_2, (-3, -195), axis=(1, 0)))*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [ ]:
 
In [356]:
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img, _ = cv2.split(img)

dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_14_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img_2, _ = cv2.split(img)

ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [357]:
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + np.roll(ch_img_2, (-3, -195), axis=(1, 0)))*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [ ]:
 
In [358]:
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_17_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img, _ = cv2.split(img)

dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_16_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img_2, _ = cv2.split(img)

ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [362]:
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + np.roll(ch_img_2, (-3, -210), axis=(1, 0)))*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [ ]:
 
In [363]:
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_19_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img, _ = cv2.split(img)

dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_18_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img_2, _ = cv2.split(img)

ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [365]:
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + np.roll(ch_img_2, (-3, -200), axis=(1, 0)))*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [ ]:
 
In [376]:
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_13_y_13_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img, _ = cv2.split(img)

dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_12_y_12_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img_2, _ = cv2.split(img)

ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [389]:
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + np.roll(ch_img_2, (177, -180), axis=(1, 0)))*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [ ]:
 
In [390]:
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img, _ = cv2.split(img)

dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_14_y_14_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img_2, _ = cv2.split(img)

ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [395]:
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + np.roll(ch_img_2, (187, -205), axis=(1, 0)))*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [ ]:
 
In [396]:
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_17_y_17_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img, _ = cv2.split(img)

dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_16_y_16_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img_2, _ = cv2.split(img)

ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + ch_img_2)*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [410]:
ch_img_1 = np.zeros_like(img)
ch_img_1[:, :, 0] = g_img

ch_img_2 = np.zeros_like(img)
ch_img_2[:, :, 1] = g_img_2

new_img = np.array((ch_img_1 + np.roll(ch_img_2, (195, -200), axis=(1, 0)))*2, dtype=np.uint8)

fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(new_img)
plt.show()
In [ ]:
 

Deductions:¶

  • Movement across x-axis is constant
  • Movement accross y-axis seems to reduce as you go farther from it
  • Movement diagonally is changing. Maybe the same as y-axis and because of it.
In [ ]:
 

Find optimal parameters for ellipses¶

In [ ]:
 
In [492]:
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = "2023_03_15_16_02_07_img_shutter_05_x_17_y_17_r_0_g_1_b_0.tiff"

img = cv2.imread(dir_path + file_path)

_, g_img, _ = cv2.split(img)
In [531]:
g_img_2 = g_img.copy()
ellipse = ((2755+390+15, 1210-390),
           (5322, 5200),
           0)

# 17,17 -> not a good fit: more elliptical
# ((2755+390+15, 1210-390),
#            (5322, 5200),
#            0)

# 13,13
# ((2755-390, 1210+390),
#            (5322, 5322),
#            0)

# 15, 17
# ((2755, 1210-380),
#            (5322, 5200),
#            0)

# 15, 13
# ((2755, 1210+400/390),
#            (5320, 5320),
#            0)

# 13,15: 
# ((2765-410, 1210-15),
#            (5320, 5320),
#            0)

# 15,15:
# ((2765, 1210),
#            (5320, 5320),
#            0)

# 17,15:
# ((2765+390, 1210-20),
#            (5320, 5320),
#            0)


cv2.ellipse(g_img_2, ellipse, 255, 4)
fig, ax = plt.subplots(1,1, figsize=(15, 15))
ax.imshow(g_img_2, cmap='gray')
plt.show()
In [ ]:
 

Observations:¶

  • For many close to the center (that are tested), a circle seems to fit the boundary of the object well
  • As the y index of the LEDs increases, it becomes tougher (impossible?) to fit a circle on the boundary of the object
In [ ]:
 
In [ ]:
 
In [8]:
def contour_points_v2(g_img):

    kernel = np.ones((15, 15), np.uint8)
    gradient = cv2.morphologyEx(g_img, cv2.MORPH_GRADIENT, kernel)
    
    blur = cv2.GaussianBlur(gradient + 20, (5, 5), 0)
    ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)

    gradient_img = np.zeros_like(g_img, dtype=np.uint8)

    contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
    if len(contours) < 1:
        return None
    
    cv2.drawContours(gradient_img, contours, -1, 255, 2)

    kernel = np.ones((25, 25),np.uint8)
    gradient_img = cv2.dilate(gradient_img, kernel, iterations = 1)    


    hull_img = np.zeros_like(g_img, dtype=np.uint8)

    blur = cv2.GaussianBlur(g_img, (3, 3), 0)
    ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)

    contours, _ = cv2.findContours(th3, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
    
    if len(contours) < 1:
        return None
    
    area_cntrs = [cv2.contourArea(cnt) for cnt in contours]
    max_cntr_index = np.argmax(area_cntrs)
    cnt = contours[max_cntr_index]

    hull = cv2.convexHull(cnt)
    cv2.drawContours(hull_img, [hull], 0, 125, 2)

    op = np.logical_and(gradient_img, hull_img).astype(np.uint8)*255

#     op[:10, :] = 0
#     op[-10:, :] = 0

#     kernel = np.ones((3, 3), np.uint8)
#     opening = cv2.morphologyEx(op, cv2.MORPH_OPEN, kernel)
    

#     final_cnt = np.take(final_cnt, list(set(np.random.randint(0, final_cnt.shape[0], size=int(final_cnt.shape[0]/10)))), axis=0)
#     imgs = np.hstack([gradient_img, hull_img, opening])
    
#     if final_cnt.shape[0] > 4:
#         return op
    
#     else:
#         return None
    return op
In [9]:
# Ivo's implemetation of Ellipse class

class Ellipse:

    C : np.array #homogeneous conic matrix x'*C*x=0

    def __init__( self, C : np.array ):

        self.C = C
        

    @classmethod
    def estimate_ellipsoid( cls, points : np.array ):
        """
        C = estimate_ellipsoid_2D ( points )
        ----------------------------------

        estimates a homogeneous ellipsoid C
        
        x'*C*x = 0
        
        from a number of point samples 
        (in a least squares sense)
        
        format of points

                [ x^1_1 . . . . .  x^N_1 ] 
                [  .                 .   ]
        points = [  .                 .   ]
                [  .                 .   ]
                [ x^1_n . . . . .  x^N_n ]
                [  1    . . . . .    1   ]

        where each point is element R^n
        and there are N points, the points  
        have to be in homogeneous coordinates.
        
        """

        n = points.shape[0]
        num_entries = int((n+1) * n/2)
        
        if points.shape[1] < num_entries-1:
            print('insufficient number of points, need at least {}'.format( num_entries-1 ) )
    
        # --------------  estimate ellipse ---------------

        # assemble linear system
        A = np.zeros( [points.shape[ 1 ], num_entries ] )

        #for all points
        for i in range( 0, points.shape[ 1 ] ):
        
            #extract point
            x = np.array([points[:, i]]).T
            #print(x)
            X  = np.matmul( x, x.T )
            #print(X)
            X2 = X.T + X - np.diag( np.diag( X ) )

            #print(X2)
            row = X2[ np.where( np.triu( np.ones( X2.shape ) ) ) ]
            #print(row)

            A[ i, : ] = row

        #print(A)

        #solve for null space 
        U,D,V = np.linalg.svd( A )

        ell = V[-1,:]


        #create conic matrix
        Cd = np.zeros( X2.shape )

        Cd[ np.where( np.triu( np.ones( X2.shape ) ) ) ] = ell

        C = Cd + Cd.T - np.diag( np.diag ( Cd ) )

        #detC =  np.linalg.det( C[0:-1,0:-1] )

        #print(detC)

        if Ellipse.is_hyperboloid( C ):
            print('Warning: result is a hyperboloid!')
        

        if Ellipse.is_paraboloid( C ):
            print('Warning: result is a paraboloid!')
        
        return C

    @classmethod
    def cov_mean_to_ellipsoid( cls, cov, mean, k ):
        """
            C = cov_mean_to_ellipsoid_homog ( Cov, Mean, k )
            ----------------------------------
        
            generates a homogeneous ellipsoid C
            
            x'*C*x = 0
        
            from covariance matrix and mean value
            of a gaussian distribution
            
            k is a confidence factor determibed by the chi^2 
            distribution. The probability that x is inside 
            the ellipsoid is P_{chi^2}(k,2). Some values 
            for k are 2.41 (70%) and k=5.99 (95%).
            
            cov  : np.array - covariance matrix (non-homogeneous)
            mean : np.array - mean value (vector)
            k : float - confidence factor, see above 

        """

        #column vector, assume we get a 1D array
        mean = np.array( [mean] ).T

        SIGMAinv = np.linalg.inv( cov )

        lastentry = np.matmul( mean.T, np.matmul( SIGMAinv, mean ) ) - k 
        C = np.vstack( [ np.hstack( [ SIGMAinv, -np.matmul(SIGMAinv, mean)] ), 
                         np.hstack( [ -np.matmul(mean.T,SIGMAinv), lastentry] ) ] )
        
        return C

    @classmethod
    def srt_to_ellipsoid( cls, S, R, T ): 
        """
         C = srt_to_ellipsoid_homog ( S, R, T )
         ----------------------------------
        
         generates a homogeneous ellipsoid C
         
         x'*C*x = 0
        
         from translation, rotation and scale parts as
        
         x'* (T^-1 * R^-1 * S^-1) ( S * R * T ) x = 0
        
         see: Ivo Ihrke: "Some Notes on Ellipses", TechReport, 2004
        """
        #function C = srt_to_ellipsoid_homog ( S, R, T )


        invT = np.linalg.inv( T )
        invR = R.T
        invS = np.linalg.inv( S )

        left  = np.matmul( invT.T, np.matmul( invR.T, invS.T ) )
        right = np.matmul( invS, np.matmul( invR, invT ) )

        C = np.matmul( left, right )
    
        return C



    @classmethod
    def decompose_ellipsoid( cls, C ):
        """
        
         [S,R,T] = decompose_ellipsoid_homog ( C )
         ----------------------------------
        
         decomposes a homogeneous ellipsoid C
         
         x'*C*x = 0
        
         into translation, rotation and scale parts
        
         x'* (T^-1 * R^-1 * S^-1) ( S * R * T ) x = 0
        
         Note that C \neq (T^-1 * R^-1 * S^-1) ( S * R * T )
         in some cases because the orientation of 
         the ellipsoid's major axis is not fixed.
         The function generates a right handed system.

        """

        #print(C)
        C_ih = C[0:-1,0:-1] 
        #print(C_ih)

        t = -np.matmul( np.linalg.pinv( C_ih ), C[0:-1,-1] )

        T = np.eye( C.shape[ 0 ] )
        T[0:-1, -1 ] = t

        U,D,V = np.linalg.svd( C_ih )

        R = np.eye( C.shape[ 0 ] )
        R[0:-1,0:-1] = V


        #S = R' * T' * C * T * R;
        S = np.matmul( R.T, np.matmul( T.T, np.matmul( C, np.matmul( T, R ) ) ) )

        S = S / S[-1,-1]

        S = np.sqrt( np.diag(1/abs(np.diag(S))))

        return S, R, T

    @classmethod
    def is_ellipsoid( cls, C ):
        """
         bool = is_ellipsoid_homog ( C )
         ----------------------------------
        
         determines whether homogeneous conic C
         
         x'*C*x = 0
         
         is an ellipsoid
        """

        #function val = conic_is_ellipsoid_homog ( C )

        detC =  np.linalg.det( C[0:-1,0:-1] )

        return detC > 0


    @classmethod
    def is_hyperboloid( cls, C ):
        """
         bool = is_hyperboloid_homog ( C )
         ----------------------------------
        
         determines whether homogeneous conic C
         
         x'*C*x = 0
         
         is a hyperboloid
        """

        
        detC =  np.linalg.det( C[0:-1,0:-1] )

        return detC < 0

    @classmethod
    def is_paraboloid( cls, C ):
        """
         bool = is_paraboloid_homog ( C )
         ----------------------------------
        
         determines whether homogeneous conic C
         
         x'*C*x = 0
         
         is a paraboloid
        """

        
        detC =  np.linalg.det( C[0:-1,0:-1] )

        #maybe with EPS ?
        return detC == 0


    def get_center( self ):
        
        S,R,T = self.decompose_ellipsoid( self.C )


        return [ T[2,0], T[2,1] ]


    def plot( self, ax, **kwargs ):
        """
        plot_ellipsoid_2D ( C, style, N )
        ----------------------------------
        
        plots a homogeneous ellipsoid C
        
        x'*C*x = 0
        
        style is a matlab plot style, e.g. 'b-' (default)
        N is the number of polygon segments to 
        approximate the ellipse
        """

        #function plot_ellipsoid_homog_2d ( C, style, N )

        if not 'style' in locals():
            style = 'b-'

        if not 'N' in locals():
            N = 100



        angles = np.linspace( -np.pi, np.pi, N )

        # generate a unit circle or a unit hyperbola
        if ( Ellipse.is_ellipsoid( self.C ) ):
            #uc  = [ cos(x); sin(x); ones(size(x)) ];
            x  = np.cos( angles )
            y  = np.sin( angles )
            w  = np.ones( angles.shape )

            uc = np.vstack([x,y,w])

        else:
            EPS = 0.1
            x   =  np.linspace( -np.pi/2+EPS, np.pi/2-EPS, int( N/2 ) )#-(pi/2-EPS):((pi-2*EPS)/(N/2)):(pi/2-EPS);
            x2  =  np.linspace( np.pi/2+EPS, 3*np.pi/2-EPS, int( N/2 ) )  #x2 =  (pi/2+EPS):((pi-2*EPS)/(N/2)):(3*pi/2-EPS);
            uc  = np.vstack( [ 1/np.cos(x), np.tan(x), np.ones( x.shape ) ] )
            uc2 = np.vstack( [ 1/np.cos(x2), np.tan(x2), np.ones( x2.shape ) ] )
        
        #decompose
        S,R,T = Ellipse.decompose_ellipsoid( self.C )

        #transform unit circle or parabola 
        transform = np.matmul( T, np.matmul( R, S ) )
        ell  =  np.matmul( transform, uc )

        if 'uc2' in locals():
            ell2  =  np.matmul( transform, uc2 )
            ax.plot( ell2[0,:], ell2[1,:], style )

   
        #plot conic
        ax.plot( ell[0,:], ell[1,:], style )

        #plot coordinate system
        tx = T[0,2]
        ty = T[1,2]
        
#         print(tx, ty)

        smin = S[0,0]
        smaj = S[1,1] 

        ax.plot( np.array( [ tx, tx+smin*R[0,0] ] ), np.array( [ ty, ty+smin*R[1,0] ] ), 'r' )
        ax.plot( np.array( [ tx, tx+smaj*R[0,1] ] ), np.array( [ ty, ty+smaj*R[1,1] ] ), 'b' )

        ax.set_aspect('equal', adjustable='box')
        return (tx, ty)
In [10]:
def fit_ellipse_3(img_name, img_8, g_img_8, ifHull=False, date_now="", method="ivo", save_path=None):
        """ Find and fit an ellipse in an image.
        it should be run on a normalized image I (normalization should exclude dead pixels which are artificially maximum)
        Parameters
        ----------
        path_dir_w_time : str
            _description_
        date_now : str
            _description_

        Returns
        -------
        Ellipse
            Ellipse in an image.
        """

        ellipsoid_points=[]
           
        img = np.array(g_img_8, dtype=np.uint8)
        op = contour_points_v2(img)
        
        op[:20, :] = 0
        op[-20:, :] = 0
        
        op[:, :20] = 0
        op[:, -20:] = 0
        
#         plt.imshow(op)
#         plt.show()
        
        x, y = np.where(op > 1)
        final_cnt = np.vstack([y, x]).T
        
        if final_cnt.shape[0] > 10000:
            final_cnt = np.take(final_cnt, list(set(np.random.randint(0, final_cnt.shape[0], size=10000))), axis=0)

        if final_cnt is None:
            return None, None, None
        else:
            ellipsoid_points = np.array(final_cnt).reshape([-1, 2])
            
        
        ellipsoid_points = np.array(ellipsoid_points).astype('float').T
        ellipsoid_points = np.vstack( [ellipsoid_points, np.ones(ellipsoid_points.shape[1])] )
#         print(ellipsoid_points)
            
        is_det_ellipse = True
        if method == "ivo":
            E = Ellipse( Ellipse.estimate_ellipsoid( ellipsoid_points ) )
            
            if Ellipse.is_ellipsoid( E.C ):
                S,R,T = Ellipse.decompose_ellipsoid(E.C)

                tx = T[0,2]
                ty = T[1,2]

                smin = S[0,0]
                smaj = S[1,1]
                
#                 f = E.C[-1, -1]
                f = np.arctan2(R[1,0], R[0, 0])
#                 print(theta)
#                 print(E.C)
            else:
                is_det_ellipse = False
        else:
            (tx, ty), (smin, smaj), f = cv2.fitEllipse(ellipsoid_points[:2, :].T.astype(np.int32))
            smin, smaj = smin/2, smaj/2


#         print((tx, ty), (smin, smaj), f)
        
        
#             
        if is_det_ellipse:
            for i in range( 0, ellipsoid_points.shape[1] ):
                draw_cross(img_8, (ellipsoid_points[0,i],ellipsoid_points[1,i]), (0, 0, 255), length=20, thickness=5)
        
            cv2.ellipse(img_8, (int(tx), int(ty)), (int(smin), int(smaj)), f, 0, 360, (255, 0, 0), 5)

            draw_cross(img_8, (tx, ty), (0, 255, 0))

#             plt.imshow(img_8)
#             plt.show()
            cv2.imshow('Resized Image', cv2.resize(img_8, (1080, 720), interpolation = cv2.INTER_LINEAR))
            cv2.waitKey(0)
            cv2.destroyAllWindows()
            
            if save_path is not None:
            
                cv2.imwrite(f'{save_path}/{img_name[:-5]}.png', img_8)
        
            return (int(tx), int(ty)), (int(smin), int(smaj)), f
        else:
            return None, None, None
In [906]:
cv2.fitEllipse(ellipsoid_points[:2, :].T.astype(np.int32))
Out[906]:
((219.3291473388672, 104.70303344726562),
 (513.478759765625, 532.2569580078125),
 53.40561294555664)
In [1054]:
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = f"2023_03_15_16_02_07_img_shutter_05_x_13_y_18_r_0_g_1_b_0.tiff"

        
img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
_, g_img, _ = cv2.split(img)
        
img_8 = cv2.convertScaleAbs(img, alpha=(255/65535))
g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
In [1055]:
center, axis, f = fit_ellipse_3(file_path, img_8, g_img_8, ifHull=True, method='ivo')
[[ 7.23449895e-07  7.32307094e-09 -1.69863270e-03]
 [ 7.32307094e-09  7.46098008e-07 -4.49513750e-04]
 [-1.69863270e-03 -4.49513750e-04 -9.99998456e-01]]
In [ ]:
 
In [11]:
def draw_cross(img, cord, color, length=50, thickness=10):
    
    x, y = int(cord[0]), int(cord[1])
    cv2.line(img, (x-length, y-length), (x+length, y+length), color, thickness)
    cv2.line(img, (x+length, y-length), (x-length, y+length), color, thickness)
    
    return img
In [ ]:
 
In [24]:
 
In [25]:
 
In [860]:
# # grid based on thumbnail images
# fit_range = np.zeros([7, 7], dtype=np.uint8)
# center_cords_ivo = np.zeros([7, 7, 2], dtype=np.float64)
# axis_lens_ivo = np.zeros([7, 7, 2], dtype=np.float64)

# for x in range(12, 19):
#     for y in range(12, 19):

#         dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
#         file_path = f"2023_03_15_16_02_07_img_shutter_05_x_{x}_y_{y}_r_0_g_1_b_0.tiff"
#         print('x:', x, 'y:', y, file_path)
        
#         img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
#         _, g_img, _ = cv2.split(img)

#         img_8 = cv2.convertScaleAbs(img, alpha=(255/65535))
#         g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
        
#         center, axis_size, f = fit_ellipse_3(file_path, img_8, g_img_8, ifHull=True)
        
#         if center is not None:
            
#             good_fit = input("Is it a good fit?")
#             fit_range[x-12, y-12] = int(good_fit)
#             center_cords_ivo[x-12, y-12] = center
#             axis_lens_ivo[x-12, y-12] = axis_size
x: 12 y: 12 2023_03_15_16_02_07_img_shutter_05_x_12_y_12_r_0_g_1_b_0.tiff
Is it a good fit?5
x: 12 y: 13 2023_03_15_16_02_07_img_shutter_05_x_12_y_13_r_0_g_1_b_0.tiff
Is it a good fit?3
x: 12 y: 14 2023_03_15_16_02_07_img_shutter_05_x_12_y_14_r_0_g_1_b_0.tiff
Is it a good fit?4
x: 12 y: 15 2023_03_15_16_02_07_img_shutter_05_x_12_y_15_r_0_g_1_b_0.tiff
Is it a good fit?2
x: 12 y: 16 2023_03_15_16_02_07_img_shutter_05_x_12_y_16_r_0_g_1_b_0.tiff
Is it a good fit?2
x: 12 y: 17 2023_03_15_16_02_07_img_shutter_05_x_12_y_17_r_0_g_1_b_0.tiff
Is it a good fit?2
x: 12 y: 18 2023_03_15_16_02_07_img_shutter_05_x_12_y_18_r_0_g_1_b_0.tiff
Warning: result is a hyperboloid!
x: 13 y: 12 2023_03_15_16_02_07_img_shutter_05_x_13_y_12_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 13 y: 13 2023_03_15_16_02_07_img_shutter_05_x_13_y_13_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 13 y: 14 2023_03_15_16_02_07_img_shutter_05_x_13_y_14_r_0_g_1_b_0.tiff
Is it a good fit?3
x: 13 y: 15 2023_03_15_16_02_07_img_shutter_05_x_13_y_15_r_0_g_1_b_0.tiff
Is it a good fit?3
x: 13 y: 16 2023_03_15_16_02_07_img_shutter_05_x_13_y_16_r_0_g_1_b_0.tiff
Warning: result is a hyperboloid!
x: 13 y: 17 2023_03_15_16_02_07_img_shutter_05_x_13_y_17_r_0_g_1_b_0.tiff
Is it a good fit?2
x: 13 y: 18 2023_03_15_16_02_07_img_shutter_05_x_13_y_18_r_0_g_1_b_0.tiff
Is it a good fit?2
x: 14 y: 12 2023_03_15_16_02_07_img_shutter_05_x_14_y_12_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 14 y: 13 2023_03_15_16_02_07_img_shutter_05_x_14_y_13_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 14 y: 14 2023_03_15_16_02_07_img_shutter_05_x_14_y_14_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 14 y: 15 2023_03_15_16_02_07_img_shutter_05_x_14_y_15_r_0_g_1_b_0.tiff
Is it a good fit?2
x: 14 y: 16 2023_03_15_16_02_07_img_shutter_05_x_14_y_16_r_0_g_1_b_0.tiff
Is it a good fit?3
x: 14 y: 17 2023_03_15_16_02_07_img_shutter_05_x_14_y_17_r_0_g_1_b_0.tiff
Is it a good fit?3
x: 14 y: 18 2023_03_15_16_02_07_img_shutter_05_x_14_y_18_r_0_g_1_b_0.tiff
Warning: result is a paraboloid!
x: 15 y: 12 2023_03_15_16_02_07_img_shutter_05_x_15_y_12_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 15 y: 13 2023_03_15_16_02_07_img_shutter_05_x_15_y_13_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 15 y: 14 2023_03_15_16_02_07_img_shutter_05_x_15_y_14_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 15 y: 15 2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff
Is it a good fit?2
x: 15 y: 16 2023_03_15_16_02_07_img_shutter_05_x_15_y_16_r_0_g_1_b_0.tiff
Is it a good fit?2
x: 15 y: 17 2023_03_15_16_02_07_img_shutter_05_x_15_y_17_r_0_g_1_b_0.tiff
Warning: result is a paraboloid!
x: 15 y: 18 2023_03_15_16_02_07_img_shutter_05_x_15_y_18_r_0_g_1_b_0.tiff
Is it a good fit?3
x: 16 y: 12 2023_03_15_16_02_07_img_shutter_05_x_16_y_12_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 16 y: 13 2023_03_15_16_02_07_img_shutter_05_x_16_y_13_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 16 y: 14 2023_03_15_16_02_07_img_shutter_05_x_16_y_14_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 16 y: 15 2023_03_15_16_02_07_img_shutter_05_x_16_y_15_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 16 y: 16 2023_03_15_16_02_07_img_shutter_05_x_16_y_16_r_0_g_1_b_0.tiff
Is it a good fit?2
x: 16 y: 17 2023_03_15_16_02_07_img_shutter_05_x_16_y_17_r_0_g_1_b_0.tiff
Is it a good fit?3
x: 16 y: 18 2023_03_15_16_02_07_img_shutter_05_x_16_y_18_r_0_g_1_b_0.tiff
Is it a good fit?3
x: 17 y: 12 2023_03_15_16_02_07_img_shutter_05_x_17_y_12_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 17 y: 13 2023_03_15_16_02_07_img_shutter_05_x_17_y_13_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 17 y: 14 2023_03_15_16_02_07_img_shutter_05_x_17_y_14_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 17 y: 15 2023_03_15_16_02_07_img_shutter_05_x_17_y_15_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 17 y: 16 2023_03_15_16_02_07_img_shutter_05_x_17_y_16_r_0_g_1_b_0.tiff
Warning: result is a paraboloid!
x: 17 y: 17 2023_03_15_16_02_07_img_shutter_05_x_17_y_17_r_0_g_1_b_0.tiff
Is it a good fit?3
x: 17 y: 18 2023_03_15_16_02_07_img_shutter_05_x_17_y_18_r_0_g_1_b_0.tiff
Is it a good fit?3
x: 18 y: 12 2023_03_15_16_02_07_img_shutter_05_x_18_y_12_r_0_g_1_b_0.tiff
Is it a good fit?5
x: 18 y: 13 2023_03_15_16_02_07_img_shutter_05_x_18_y_13_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 18 y: 14 2023_03_15_16_02_07_img_shutter_05_x_18_y_14_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 18 y: 15 2023_03_15_16_02_07_img_shutter_05_x_18_y_15_r_0_g_1_b_0.tiff
Warning: result is a paraboloid!
x: 18 y: 16 2023_03_15_16_02_07_img_shutter_05_x_18_y_16_r_0_g_1_b_0.tiff
Is it a good fit?3
x: 18 y: 17 2023_03_15_16_02_07_img_shutter_05_x_18_y_17_r_0_g_1_b_0.tiff
Is it a good fit?2
x: 18 y: 18 2023_03_15_16_02_07_img_shutter_05_x_18_y_18_r_0_g_1_b_0.tiff
Is it a good fit?5
In [890]:
plt.figure(figsize=(15, 10))

center_leds = center_cords_ivo[3, 3]
y_unit_len = center_cords_ivo[3, 2][1] - center_cords_ivo[3, 3][1] 
x_unit_len = center_cords_ivo[4, 3][0] - center_cords_ivo[3, 3][0]

center_leds, y_unit_len, x_unit_len

x_axis = [center_leds[0] - i*x_unit_len for i in range(-6, 6)]
y_axis = [center_leds[1] - i*y_unit_len for i in range(-5, 5)]

xx, yy = np.meshgrid(x_axis, y_axis)
plt.scatter(xx, yy, marker=".", color='black')

x_cords = [x if x != 0.0 else np.nan for x in center_cords_ivo.reshape([-1, 2])[:, 0]]
y_cords = [y if y != 0.0 else np.nan for y in center_cords_ivo.reshape([-1, 2])[:, 1]]
fit_range_2 = [y if y != 0 else np.nan for y in fit_range.reshape(-1)]


plt.scatter(x_cords, y_cords, marker="*", c=fit_range_2)
# plt.grid()

count = 0
for x in range(12, 19):
    for y in range(12, 19):
        plt.annotate(f"({x}, {y})",(x_cords[count], y_cords[count]), fontsize=10)
        count += 1

plt.title("Center coordinates of Ellipse Fit (Ivo) and LED Position Annotations")
# plt.xlim(210, 350)
# plt.ylim(40, 200)
plt.colorbar()
plt.show()
In [877]:
 
In [884]:
 
In [894]:
axis_lens_pro = axis_lens_ivo.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan

fig, ax = plt.subplots(1, 2, figsize=(18, 8), layout='tight')

xx, yy = np.meshgrid([i for i in range(12, 19)], [i for i in range(12, 19)])

cp = ax[0].imshow(axis_lens_pro[:, :, 0].T, interpolation='none',
                  extent=[12, 19, 12, 19]
                 )
fig.colorbar(cp)
ax[0].grid()
ax[0].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[0].set_title("Minor Axis Lengths")
ax[0].set_xlabel('LED Positions')
ax[0].set_ylabel('LED Positions')


cp = ax[1].imshow(axis_lens_pro[:, :, 1].T, interpolation='none', 
                 extent=[12, 19, 12, 19]
                 )
fig.colorbar(cp)
ax[1].grid()
ax[1].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[1].set_title("Major Axis Lengths")
ax[1].set_xlabel('LED Positions')
ax[1].set_ylabel('LED Positions')

plt.show()
In [1178]:
axis_lens_pro = axis_lens_ivo.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan

plt.hist(axis_lens_pro[:, :, 0].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.hist(axis_lens_pro[:, :, 1].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.legend()
plt.show()
In [ ]:
 
In [970]:
# # grid based on thumbnail images
# fit_range_cv2 = np.zeros([7, 7], dtype=np.uint8)
# center_cords_cv2 = np.zeros([7, 7, 2], dtype=np.float64)
# axis_lens_cv2 = np.zeros([7, 7, 2], dtype=np.float64)

# save_dir = "C:/Users/Kazim/Desktop/Overall/VignettingCircleDetection/310823_ellipse_cv2_fit_tiff_size"
# for x in range(12, 19):
#     for y in range(12, 19):

#         dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
#         file_path = f"2023_03_15_16_02_07_img_shutter_05_x_{x}_y_{y}_r_0_g_1_b_0.tiff"
#         print('x:', x, 'y:', y, file_path)
        
#         img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
#         _, g_img, _ = cv2.split(img)

#         img_8 = cv2.convertScaleAbs(img, alpha=(255/65535))
#         g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
        
#         center, axis_size, f = fit_ellipse_3(file_path, img_8, g_img_8, ifHull=True, method='opencv', save_path=save_dir)
        
#         if center is not None:
            
#             good_fit = input("Is it a good fit?")
#             fit_range_cv2[x-12, y-12] = int(good_fit)
#             center_cords_cv2[x-12, y-12] = center
#             axis_lens_cv2[x-12, y-12] = axis_size
x: 12 y: 12 2023_03_15_16_02_07_img_shutter_05_x_12_y_12_r_0_g_1_b_0.tiff
Is it a good fit?5
x: 12 y: 13 2023_03_15_16_02_07_img_shutter_05_x_12_y_13_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 12 y: 14 2023_03_15_16_02_07_img_shutter_05_x_12_y_14_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 12 y: 15 2023_03_15_16_02_07_img_shutter_05_x_12_y_15_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 12 y: 16 2023_03_15_16_02_07_img_shutter_05_x_12_y_16_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 12 y: 17 2023_03_15_16_02_07_img_shutter_05_x_12_y_17_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 12 y: 18 2023_03_15_16_02_07_img_shutter_05_x_12_y_18_r_0_g_1_b_0.tiff
Is it a good fit?2
x: 13 y: 12 2023_03_15_16_02_07_img_shutter_05_x_13_y_12_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 13 y: 13 2023_03_15_16_02_07_img_shutter_05_x_13_y_13_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 13 y: 14 2023_03_15_16_02_07_img_shutter_05_x_13_y_14_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 13 y: 15 2023_03_15_16_02_07_img_shutter_05_x_13_y_15_r_0_g_1_b_0.tiff
Is it a good fit?2
x: 13 y: 16 2023_03_15_16_02_07_img_shutter_05_x_13_y_16_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 13 y: 17 2023_03_15_16_02_07_img_shutter_05_x_13_y_17_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 13 y: 18 2023_03_15_16_02_07_img_shutter_05_x_13_y_18_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 14 y: 12 2023_03_15_16_02_07_img_shutter_05_x_14_y_12_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 14 y: 13 2023_03_15_16_02_07_img_shutter_05_x_14_y_13_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 14 y: 14 2023_03_15_16_02_07_img_shutter_05_x_14_y_14_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 14 y: 15 2023_03_15_16_02_07_img_shutter_05_x_14_y_15_r_0_g_1_b_0.tiff
Is it a good fit?2
x: 14 y: 16 2023_03_15_16_02_07_img_shutter_05_x_14_y_16_r_0_g_1_b_0.tiff
Is it a good fit?2
x: 14 y: 17 2023_03_15_16_02_07_img_shutter_05_x_14_y_17_r_0_g_1_b_0.tiff
Is it a good fit?2
x: 14 y: 18 2023_03_15_16_02_07_img_shutter_05_x_14_y_18_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 15 y: 12 2023_03_15_16_02_07_img_shutter_05_x_15_y_12_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 15 y: 13 2023_03_15_16_02_07_img_shutter_05_x_15_y_13_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 15 y: 14 2023_03_15_16_02_07_img_shutter_05_x_15_y_14_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 15 y: 15 2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 15 y: 16 2023_03_15_16_02_07_img_shutter_05_x_15_y_16_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 15 y: 17 2023_03_15_16_02_07_img_shutter_05_x_15_y_17_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 15 y: 18 2023_03_15_16_02_07_img_shutter_05_x_15_y_18_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 16 y: 12 2023_03_15_16_02_07_img_shutter_05_x_16_y_12_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 16 y: 13 2023_03_15_16_02_07_img_shutter_05_x_16_y_13_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 16 y: 14 2023_03_15_16_02_07_img_shutter_05_x_16_y_14_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 16 y: 15 2023_03_15_16_02_07_img_shutter_05_x_16_y_15_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 16 y: 16 2023_03_15_16_02_07_img_shutter_05_x_16_y_16_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 16 y: 17 2023_03_15_16_02_07_img_shutter_05_x_16_y_17_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 16 y: 18 2023_03_15_16_02_07_img_shutter_05_x_16_y_18_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 17 y: 12 2023_03_15_16_02_07_img_shutter_05_x_17_y_12_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 17 y: 13 2023_03_15_16_02_07_img_shutter_05_x_17_y_13_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 17 y: 14 2023_03_15_16_02_07_img_shutter_05_x_17_y_14_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 17 y: 15 2023_03_15_16_02_07_img_shutter_05_x_17_y_15_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 17 y: 16 2023_03_15_16_02_07_img_shutter_05_x_17_y_16_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 17 y: 17 2023_03_15_16_02_07_img_shutter_05_x_17_y_17_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 17 y: 18 2023_03_15_16_02_07_img_shutter_05_x_17_y_18_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 18 y: 12 2023_03_15_16_02_07_img_shutter_05_x_18_y_12_r_0_g_1_b_0.tiff
Is it a good fit?5
x: 18 y: 13 2023_03_15_16_02_07_img_shutter_05_x_18_y_13_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 18 y: 14 2023_03_15_16_02_07_img_shutter_05_x_18_y_14_r_0_g_1_b_0.tiff
Is it a good fit?2
x: 18 y: 15 2023_03_15_16_02_07_img_shutter_05_x_18_y_15_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 18 y: 16 2023_03_15_16_02_07_img_shutter_05_x_18_y_16_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 18 y: 17 2023_03_15_16_02_07_img_shutter_05_x_18_y_17_r_0_g_1_b_0.tiff
Is it a good fit?1
x: 18 y: 18 2023_03_15_16_02_07_img_shutter_05_x_18_y_18_r_0_g_1_b_0.tiff
Is it a good fit?5
In [972]:
plt.figure(figsize=(15, 10))

center_leds = center_cords_cv2[3, 3]
y_unit_len = center_cords_cv2[3, 2][1] - center_cords_cv2[3, 3][1] 
x_unit_len = center_cords_cv2[4, 3][0] - center_cords_cv2[3, 3][0]

center_leds, y_unit_len, x_unit_len

x_axis = [center_leds[0] - i*x_unit_len for i in range(-4, 4)]
y_axis = [center_leds[1] - i*y_unit_len for i in range(-4, 4)]

xx, yy = np.meshgrid(x_axis, y_axis)
plt.scatter(xx, yy, marker=".", color='black')

x_cords = [x if x != 0.0 else np.nan for x in center_cords_cv2.reshape([-1, 2])[:, 0]]
y_cords = [y if y != 0.0 else np.nan for y in center_cords_cv2.reshape([-1, 2])[:, 1]]
fit_range_2 = [y if y != 0 else np.nan for y in fit_range_cv2.reshape(-1)]


plt.scatter(x_cords, y_cords, marker="*", c=fit_range_2)
# plt.grid()

count = 0
for x in range(12, 19):
    for y in range(12, 19):
        plt.annotate(f"({x}, {y})",(x_cords[count], y_cords[count]), fontsize=10)
        count += 1

plt.title("Center coordinates of Ellipse Fit (Ivo) and LED Position Annotations")
# plt.xlim(210, 350)
# plt.ylim(40, 200)
plt.colorbar()
plt.show()
In [973]:
axis_lens_pro = axis_lens_cv2.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan

fig, ax = plt.subplots(1, 2, figsize=(18, 8), layout='tight')

xx, yy = np.meshgrid([i for i in range(12, 19)], [i for i in range(12, 19)])

cp = ax[0].imshow(axis_lens_pro[:, :, 0].T, interpolation='none',
                  extent=[12, 19, 12, 19]
                 )
fig.colorbar(cp)
ax[0].grid()
ax[0].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[0].set_title("Minor Axis Lengths")
ax[0].set_xlabel('LED Positions')
ax[0].set_ylabel('LED Positions')


cp = ax[1].imshow(axis_lens_pro[:, :, 1].T, interpolation='none', 
                 extent=[12, 19, 12, 19]
                 )
fig.colorbar(cp)
ax[1].grid()
ax[1].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[1].set_title("Major Axis Lengths")
ax[1].set_xlabel('LED Positions')
ax[1].set_ylabel('LED Positions')

plt.show()
In [ ]:
 
In [1179]:
axis_lens_pro = axis_lens_cv2.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan

plt.hist(axis_lens_pro[:, :, 0].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.hist(axis_lens_pro[:, :, 1].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.legend()
plt.show()
In [ ]:
# fixing A, B, C, F
In [977]:
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = f"2023_03_15_16_02_07_img_shutter_05_x_16_y_15_r_0_g_1_b_0.tiff"

        
img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
_, g_img, _ = cv2.split(img)
        
img_8 = cv2.convertScaleAbs(img, alpha=(255/65535))
g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
In [978]:
center, axis, f = fit_ellipse_3(file_path, img_8, g_img_8, ifHull=True, method='ivo') # good fit
[[-3.24326055e-07  4.75306325e-10  9.59646019e-04]
 [ 4.75306325e-10 -3.27813554e-07  3.93645616e-04]
 [ 9.59646019e-04  3.93645616e-04 -9.99999462e-01]]
In [ ]:
np.cos()
In [979]:
img = np.array(g_img_8, dtype=np.uint8)
op = contour_points_v2(img)
        
op[:20, :] = 0
op[-20:, :] = 0
        
op[:, :20] = 0
op[:, -20:] = 0

        
x, y = np.where(op > 1)
final_cnt = np.vstack([y, x]).T
        
if final_cnt.shape[0] > 10000:
    final_cnt = np.take(final_cnt, list(set(np.random.randint(0, final_cnt.shape[0], size=10000))), axis=0)

ellipsoid_points = np.array(final_cnt).reshape([-1, 2])
ellipsoid_points = np.array(ellipsoid_points).astype('float').T
ellipsoid_points
Out[979]:
array([[ 529.,  530.,  531., ..., 2178., 2179., 2180.],
       [  93.,   93.,   93., ..., 3732., 3732., 3732.]])
In [987]:
b = -3.24326055e-07*x**2 + -3.27813554e-07*y**2 + f
A = np.array([x, y]).T
A, -b
Out[987]:
(array([[  93,  529],
        [  93,  530],
        [  93,  531],
        ...,
        [3732, 2178],
        [3732, 2179],
        [3732, 2180]], dtype=int64),
 array([1.09454023, 1.09488739, 1.0952352 , ..., 7.07219959, 7.07362787,
        7.07505681]))
In [1012]:
np.linalg.lstsq(A, -b, rcond=None)
Out[1012]:
(array([0.00077151, 0.00194168]),
 array([0.65526845]),
 2,
 array([338927.99799414, 121566.05420793]))
In [991]:
9.59646019e-04*2, 3.93645616e-04*2
Out[991]:
(0.001919292038, 0.000787291232)
In [1011]:
9.26610761e-04*2, -9.39963645e-04*2
Out[1011]:
(0.001853221522, -0.00187992729)
In [1072]:
b = 4.75306325e-10*x*y + -9.99999462e-01
A = np.array([x**2, y**2, x, y]).T
A, -b
Out[1072]:
(array([[    8649,   279841,       93,      529],
        [    8649,   280900,       93,      530],
        [    8649,   281961,       93,      531],
        ...,
        [13927824,  4743684,     3732,     2178],
        [13927824,  4748041,     3732,     2179],
        [13927824,  4752400,     3732,     2180]], dtype=int64),
 array([0.99997608, 0.99997603, 0.99997599, ..., 0.99613603, 0.99613426,
        0.99613248]))
In [1085]:
A, B, D, E = np.linalg.lstsq(A, -b, rcond=None)[0]
In [1087]:
A, B, D/2, E/2
Out[1087]:
(-3.268969844369404e-07,
 -3.236289044146954e-07,
 0.0003930310412070922,
 0.0009580904040238742)
In [ ]:
 
In [ ]:
 
In [12]:
def fit_ellipse_4(img_name, img_8, g_img_8, save_path=None):
        """ Find ellipse center position vector where A=-3.24326055e-07, B=0, C=-3.27813554e-07, F=-9.99999462e-01
        Parameters
        ----------
        Returns
        -------
        Ellipse
            Ellipse in an image.
        """

        ellipsoid_points=[]
           
        img = np.array(g_img_8, dtype=np.uint8)
        op = contour_points_v2(img)
        
        op[:20, :] = 0
        op[-20:, :] = 0
        
        op[:, :20] = 0
        op[:, -20:] = 0

        x, y = np.where(op > 0)
        final_cnt = np.vstack([y, x]).T

        if final_cnt is None:
            return None, None, None
        else:
            if final_cnt.shape[0] > 10000:
                final_cnt = np.take(final_cnt, list(set(np.random.randint(0, final_cnt.shape[0], size=10000))), axis=0)
            ellipsoid_points = np.array(final_cnt).reshape([-1, 2])
            
        
        ellipsoid_points = np.array(ellipsoid_points).astype('float').T
        ellipsoid_points = np.vstack( [ellipsoid_points, np.ones(ellipsoid_points.shape[1])] )
            
        is_det_ellipse = True
        
#         C = -3.24326055e-07
        B = 4.75306325e-10
#         A = -3.27813554e-07
        F = -9.99999462e-01*2
#         b = A*x**2 + B*x*y + C*y**2 + F
        b = B*x*y + F

#         eqs = np.array([x**2, y**2, x, y]).T
        
#         C, A, E, D = np.linalg.lstsq(eqs, -b, rcond=None)[0]/2
        
#         matC = np.array([[A, B, D/2], [B, C, E/2], [D/2, E/2, F]])

        eqs = np.array([x**2 + y**2, x, y]).T
        
        A, E, D = np.linalg.lstsq(eqs, -b, rcond=None)[0]/2
        
        matC = np.array([[A, B, D/2], [B, A, E/2], [D/2, E/2, F/2]])

        E = Ellipse( matC )
#         print(Ellipse.is_ellipsoid( E.C ))
        print(matC)
#         print(np.linalg.lstsq(eqs, -b, rcond=None))
            
        if Ellipse.is_ellipsoid( E.C ):
            S,R,T = Ellipse.decompose_ellipsoid(E.C)

            tx = T[0,2]
            ty = T[1,2]

            smin = S[0,0]
            smaj = S[1,1]
            
            theta = np.arctan2(R[1,0], R[0, 0])
#             print(theta)
            
                
#             f = E.C[-1, -1]
#             print(E.C)
        else:
            is_det_ellipse = False
             
        if is_det_ellipse:
            for i in range( 0, ellipsoid_points.shape[1] ):
                draw_cross(img_8, (ellipsoid_points[0,i],ellipsoid_points[1,i]), (0, 0, 255), length=20, thickness=5)
        
            cv2.ellipse(img_8, (int(tx), int(ty)), (int(smin), int(smaj)), theta, 0, 360, (255, 0, 0), 5)

            draw_cross(img_8, (tx, ty), (0, 255, 0))

#             plt.imshow(img_8)
#             plt.show()
            cv2.imshow('Resized Image', cv2.resize(img_8, (1080, 720), interpolation = cv2.INTER_LINEAR))
            cv2.waitKey(0)
            cv2.destroyAllWindows()
            
            if save_path is not None:
            
                cv2.imwrite(f'{save_path}/{img_name[:-5]}.png', img_8)
        
            return (int(tx), int(ty)), (int(smin), int(smaj)), theta
        else:
            return None, None, None
In [1166]:
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = f"2023_03_15_16_02_07_img_shutter_05_x_15_y_17_r_0_g_1_b_0.tiff"

        
img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
_, g_img, _ = cv2.split(img)
        
img_8 = cv2.convertScaleAbs(img, alpha=(255/65535))
g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
In [1167]:
center, axis_size, f = fit_ellipse_4(file_path, img_8, g_img_8, save_path=None)
[[-9.96262090e-07  4.75306325e-10  2.74601838e-03]
 [ 4.75306325e-10 -9.96262090e-07  7.62551767e-04]
 [ 2.74601838e-03  7.62551767e-04 -9.99999462e-01]]
In [ ]:
 
In [1169]:
# # grid based on thumbnail images
# fit_range_circle = np.zeros([7, 7], dtype=np.uint8)
# center_cords_circle = np.zeros([7, 7, 2], dtype=np.float64)
# axis_lens_circle = np.zeros([7, 7, 2], dtype=np.float64)

# save_dir = "C:/Users/Kazim/Desktop/Overall/VignettingCircleDetection/310823_circle_fit_tiff_size"
# for x in range(12, 19):
#     for y in range(12, 19):

#         dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
#         file_path = f"2023_03_15_16_02_07_img_shutter_05_x_{x}_y_{y}_r_0_g_1_b_0.tiff"
#         print('x:', x, 'y:', y, file_path)
        
#         img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
#         _, g_img, _ = cv2.split(img)

#         img_8 = cv2.convertScaleAbs(img, alpha=(255/65535))
#         g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
        
#         center, axis_size, f = fit_ellipse_4(file_path, img_8, g_img_8, save_path=save_dir)
        
#         if center is not None:
            
#             good_fit = input("Is it a good fit?")
#             fit_range_circle[x-12, y-12] = int(good_fit)
#             center_cords_circle[x-12, y-12] = center
#             axis_lens_circle[x-12, y-12] = axis_size
x: 12 y: 12 2023_03_15_16_02_07_img_shutter_05_x_12_y_12_r_0_g_1_b_0.tiff
[[-1.51944827e-07  4.75306325e-10  4.28640374e-04]
 [ 4.75306325e-10 -1.51944827e-07  2.64667160e-04]
 [ 4.28640374e-04  2.64667160e-04 -9.99999462e-01]]
Is it a good fit?5
x: 12 y: 13 2023_03_15_16_02_07_img_shutter_05_x_12_y_13_r_0_g_1_b_0.tiff
[[-8.05421194e-07  4.75306325e-10  1.84482464e-03]
 [ 4.75306325e-10 -8.05421194e-07  1.29403723e-03]
 [ 1.84482464e-03  1.29403723e-03 -9.99999462e-01]]
Is it a good fit?3
x: 12 y: 14 2023_03_15_16_02_07_img_shutter_05_x_12_y_14_r_0_g_1_b_0.tiff
[[ 1.82876905e-06  4.75306325e-10 -3.95544562e-03]
 [ 4.75306325e-10  1.82876905e-06 -2.53086219e-03]
 [-3.95544562e-03 -2.53086219e-03 -9.99999462e-01]]
Is it a good fit?1
x: 12 y: 15 2023_03_15_16_02_07_img_shutter_05_x_12_y_15_r_0_g_1_b_0.tiff
[[ 8.30978298e-07  4.75306325e-10 -1.78516376e-03]
 [ 4.75306325e-10  8.30978298e-07 -9.72636790e-04]
 [-1.78516376e-03 -9.72636790e-04 -9.99999462e-01]]
Is it a good fit?2
x: 12 y: 16 2023_03_15_16_02_07_img_shutter_05_x_12_y_16_r_0_g_1_b_0.tiff
[[ 5.53879588e-07  4.75306325e-10 -1.18126017e-03]
 [ 4.75306325e-10  5.53879588e-07 -5.30531858e-04]
 [-1.18126017e-03 -5.30531858e-04 -9.99999462e-01]]
Is it a good fit?2
x: 12 y: 17 2023_03_15_16_02_07_img_shutter_05_x_12_y_17_r_0_g_1_b_0.tiff
[[ 4.35008066e-07  4.75306325e-10 -9.22324007e-04]
 [ 4.75306325e-10  4.35008066e-07 -3.19552659e-04]
 [-9.22324007e-04 -3.19552659e-04 -9.99999462e-01]]
Is it a good fit?2
x: 12 y: 18 2023_03_15_16_02_07_img_shutter_05_x_12_y_18_r_0_g_1_b_0.tiff
[[ 3.82984548e-07  4.75306325e-10 -8.12447687e-04]
 [ 4.75306325e-10  3.82984548e-07 -1.98597659e-04]
 [-8.12447687e-04 -1.98597659e-04 -9.99999462e-01]]
Is it a good fit?2
x: 13 y: 12 2023_03_15_16_02_07_img_shutter_05_x_13_y_12_r_0_g_1_b_0.tiff
[[-5.58856538e-07  4.75306325e-10  1.33007273e-03]
 [ 4.75306325e-10 -5.58856538e-07  1.00255279e-03]
 [ 1.33007273e-03  1.00255279e-03 -9.99999462e-01]]
Is it a good fit?1
x: 13 y: 13 2023_03_15_16_02_07_img_shutter_05_x_13_y_13_r_0_g_1_b_0.tiff
[[-9.23869595e-07  4.75306325e-10  2.18874479e-03]
 [ 4.75306325e-10 -9.23869595e-07  1.48275359e-03]
 [ 2.18874479e-03  1.48275359e-03 -9.99999462e-01]]
Is it a good fit?1
x: 13 y: 14 2023_03_15_16_02_07_img_shutter_05_x_13_y_14_r_0_g_1_b_0.tiff
[[-2.27036215e-06  4.75306325e-10  5.37034775e-03]
 [ 4.75306325e-10 -2.27036215e-06  3.18468554e-03]
 [ 5.37034775e-03  3.18468554e-03 -9.99999462e-01]]
Is it a good fit?1
x: 13 y: 15 2023_03_15_16_02_07_img_shutter_05_x_13_y_15_r_0_g_1_b_0.tiff
[[ 3.06726902e-06  4.75306325e-10 -7.19762759e-03]
 [ 4.75306325e-10  3.06726902e-06 -3.60634197e-03]
 [-7.19762759e-03 -3.60634197e-03 -9.99999462e-01]]
Is it a good fit?2
x: 13 y: 16 2023_03_15_16_02_07_img_shutter_05_x_13_y_16_r_0_g_1_b_0.tiff
[[ 1.10581527e-06  4.75306325e-10 -2.58177670e-03]
 [ 4.75306325e-10  1.10581527e-06 -1.05802929e-03]
 [-2.58177670e-03 -1.05802929e-03 -9.99999462e-01]]
Is it a good fit?2
x: 13 y: 17 2023_03_15_16_02_07_img_shutter_05_x_13_y_17_r_0_g_1_b_0.tiff
[[ 7.32686924e-07  4.75306325e-10 -1.70581000e-03]
 [ 4.75306325e-10  7.32686924e-07 -5.39313733e-04]
 [-1.70581000e-03 -5.39313733e-04 -9.99999462e-01]]
Is it a good fit?2
x: 13 y: 18 2023_03_15_16_02_07_img_shutter_05_x_13_y_18_r_0_g_1_b_0.tiff
[[ 5.85509650e-07  4.75306325e-10 -1.36058178e-03]
 [ 4.75306325e-10  5.85509650e-07 -3.02023162e-04]
 [-1.36058178e-03 -3.02023162e-04 -9.99999462e-01]]
Is it a good fit?2
x: 14 y: 12 2023_03_15_16_02_07_img_shutter_05_x_14_y_12_r_0_g_1_b_0.tiff
[[-3.66265272e-07  4.75306325e-10  9.41946218e-04]
 [ 4.75306325e-10 -3.66265272e-07  6.56238900e-04]
 [ 9.41946218e-04  6.56238900e-04 -9.99999462e-01]]
Is it a good fit?1
x: 14 y: 13 2023_03_15_16_02_07_img_shutter_05_x_14_y_13_r_0_g_1_b_0.tiff
[[-4.88584245e-07  4.75306325e-10  1.25433906e-03]
 [ 4.75306325e-10 -4.88584245e-07  7.83008428e-04]
 [ 1.25433906e-03  7.83008428e-04 -9.99999462e-01]]
Is it a good fit?1
x: 14 y: 14 2023_03_15_16_02_07_img_shutter_05_x_14_y_14_r_0_g_1_b_0.tiff
[[-7.02470516e-07  4.75306325e-10  1.80119649e-03]
 [ 4.75306325e-10 -7.02470516e-07  9.88621422e-04]
 [ 1.80119649e-03  9.88621422e-04 -9.99999462e-01]]
Is it a good fit?1
x: 14 y: 15 2023_03_15_16_02_07_img_shutter_05_x_14_y_15_r_0_g_1_b_0.tiff
[[-1.26810745e-06  4.75306325e-10  3.23967294e-03]
 [ 4.75306325e-10 -1.26810745e-06  1.51545546e-03]
 [ 3.23967294e-03  1.51545546e-03 -9.99999462e-01]]
Is it a good fit?2
x: 14 y: 16 2023_03_15_16_02_07_img_shutter_05_x_14_y_16_r_0_g_1_b_0.tiff
[[-3.00324635e-06  4.75306325e-10  7.67149042e-03]
 [ 4.75306325e-10 -3.00324635e-06  2.95940159e-03]
 [ 7.67149042e-03  2.95940159e-03 -9.99999462e-01]]
Is it a good fit?2
x: 14 y: 17 2023_03_15_16_02_07_img_shutter_05_x_14_y_17_r_0_g_1_b_0.tiff
[[ 2.78928085e-06  4.75306325e-10 -7.07337439e-03]
 [ 4.75306325e-10  2.78928085e-06 -2.04190106e-03]
 [-7.07337439e-03 -2.04190106e-03 -9.99999462e-01]]
Is it a good fit?2
x: 14 y: 18 2023_03_15_16_02_07_img_shutter_05_x_14_y_18_r_0_g_1_b_0.tiff
[[ 1.60175355e-06  4.75306325e-10 -4.06730316e-03]
 [ 4.75306325e-10  1.60175355e-06 -8.32599901e-04]
 [-4.06730316e-03 -8.32599901e-04 -9.99999462e-01]]
Is it a good fit?2
x: 15 y: 12 2023_03_15_16_02_07_img_shutter_05_x_15_y_12_r_0_g_1_b_0.tiff
[[-2.66292824e-07  4.75306325e-10  7.36336097e-04]
 [ 4.75306325e-10 -2.66292824e-07  4.76573016e-04]
 [ 7.36336097e-04  4.76573016e-04 -9.99999462e-01]]
Is it a good fit?1
x: 15 y: 13 2023_03_15_16_02_07_img_shutter_05_x_15_y_13_r_0_g_1_b_0.tiff
[[-3.25284852e-07  4.75306325e-10  8.98235317e-04]
 [ 4.75306325e-10 -3.25284852e-07  5.20749145e-04]
 [ 8.98235317e-04  5.20749145e-04 -9.99999462e-01]]
Is it a good fit?1
x: 15 y: 14 2023_03_15_16_02_07_img_shutter_05_x_15_y_14_r_0_g_1_b_0.tiff
[[-4.03746681e-07  4.75306325e-10  1.11518083e-03]
 [ 4.75306325e-10 -4.03746681e-07  5.67327951e-04]
 [ 1.11518083e-03  5.67327951e-04 -9.99999462e-01]]
Is it a good fit?1
x: 15 y: 15 2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff
[[-5.21047365e-07  4.75306325e-10  1.43702280e-03]
 [ 4.75306325e-10 -5.21047365e-07  6.26811405e-04]
 [ 1.43702280e-03  6.26811405e-04 -9.99999462e-01]]
Is it a good fit?1
x: 15 y: 16 2023_03_15_16_02_07_img_shutter_05_x_15_y_16_r_0_g_1_b_0.tiff
[[-7.03250637e-07  4.75306325e-10  1.93891868e-03]
 [ 4.75306325e-10 -7.03250637e-07  6.94528824e-04]
 [ 1.93891868e-03  6.94528824e-04 -9.99999462e-01]]
Is it a good fit?2
x: 15 y: 17 2023_03_15_16_02_07_img_shutter_05_x_15_y_17_r_0_g_1_b_0.tiff
[[-9.96262090e-07  4.75306325e-10  2.74601838e-03]
 [ 4.75306325e-10 -9.96262090e-07  7.62551767e-04]
 [ 2.74601838e-03  7.62551767e-04 -9.99999462e-01]]
Is it a good fit?2
x: 15 y: 18 2023_03_15_16_02_07_img_shutter_05_x_15_y_18_r_0_g_1_b_0.tiff
[[-1.36544789e-06  4.75306325e-10  3.75679371e-03]
 [ 4.75306325e-10 -1.36544789e-06  7.67007660e-04]
 [ 3.75679371e-03  7.67007660e-04 -9.99999462e-01]]
Is it a good fit?2
x: 16 y: 12 2023_03_15_16_02_07_img_shutter_05_x_16_y_12_r_0_g_1_b_0.tiff
[[-2.06004470e-07  4.75306325e-10  6.09752267e-04]
 [ 4.75306325e-10 -2.06004470e-07  3.67305162e-04]
 [ 6.09752267e-04  3.67305162e-04 -9.99999462e-01]]
Is it a good fit?1
x: 16 y: 13 2023_03_15_16_02_07_img_shutter_05_x_16_y_13_r_0_g_1_b_0.tiff
[[-2.38580862e-07  4.75306325e-10  7.06047598e-04]
 [ 4.75306325e-10 -2.38580862e-07  3.80434787e-04]
 [ 7.06047598e-04  3.80434787e-04 -9.99999462e-01]]
Is it a good fit?2
x: 16 y: 14 2023_03_15_16_02_07_img_shutter_05_x_16_y_14_r_0_g_1_b_0.tiff
[[-2.77778132e-07  4.75306325e-10  8.22098410e-04]
 [ 4.75306325e-10 -2.77778132e-07  3.89069642e-04]
 [ 8.22098410e-04  3.89069642e-04 -9.99999462e-01]]
Is it a good fit?1
x: 16 y: 15 2023_03_15_16_02_07_img_shutter_05_x_16_y_15_r_0_g_1_b_0.tiff
[[-3.25928641e-07  4.75306325e-10  9.65053368e-04]
 [ 4.75306325e-10 -3.25928641e-07  3.89726595e-04]
 [ 9.65053368e-04  3.89726595e-04 -9.99999462e-01]]
Is it a good fit?1
x: 16 y: 16 2023_03_15_16_02_07_img_shutter_05_x_16_y_16_r_0_g_1_b_0.tiff
[[-3.89253374e-07  4.75306325e-10  1.15176103e-03]
 [ 4.75306325e-10 -3.89253374e-07  3.81414201e-04]
 [ 1.15176103e-03  3.81414201e-04 -9.99999462e-01]]
Is it a good fit?2
x: 16 y: 17 2023_03_15_16_02_07_img_shutter_05_x_16_y_17_r_0_g_1_b_0.tiff
[[-4.61463753e-07  4.75306325e-10  1.36490746e-03]
 [ 4.75306325e-10 -4.61463753e-07  3.52580456e-04]
 [ 1.36490746e-03  3.52580456e-04 -9.99999462e-01]]
Is it a good fit?2
x: 16 y: 18 2023_03_15_16_02_07_img_shutter_05_x_16_y_18_r_0_g_1_b_0.tiff
[[-5.47788559e-07  4.75306325e-10  1.61837109e-03]
 [ 4.75306325e-10 -5.47788559e-07  2.99621247e-04]
 [ 1.61837109e-03  2.99621247e-04 -9.99999462e-01]]
Is it a good fit?3
x: 17 y: 12 2023_03_15_16_02_07_img_shutter_05_x_17_y_12_r_0_g_1_b_0.tiff
[[-1.64699627e-07  4.75306325e-10  5.20141144e-04]
 [ 4.75306325e-10 -1.64699627e-07  2.94129674e-04]
 [ 5.20141144e-04  2.94129674e-04 -9.99999462e-01]]
Is it a good fit?2
x: 17 y: 13 2023_03_15_16_02_07_img_shutter_05_x_17_y_13_r_0_g_1_b_0.tiff
[[-1.84993418e-07  4.75306325e-10  5.83879266e-04]
 [ 4.75306325e-10 -1.84993418e-07  2.95225429e-04]
 [ 5.83879266e-04  2.95225429e-04 -9.99999462e-01]]
Is it a good fit?2
x: 17 y: 14 2023_03_15_16_02_07_img_shutter_05_x_17_y_14_r_0_g_1_b_0.tiff
[[-2.07921155e-07  4.75306325e-10  6.56266515e-04]
 [ 4.75306325e-10 -2.07921155e-07  2.90163246e-04]
 [ 6.56266515e-04  2.90163246e-04 -9.99999462e-01]]
Is it a good fit?2
x: 17 y: 15 2023_03_15_16_02_07_img_shutter_05_x_17_y_15_r_0_g_1_b_0.tiff
[[-2.34890378e-07  4.75306325e-10  7.42205549e-04]
 [ 4.75306325e-10 -2.34890378e-07  2.78283250e-04]
 [ 7.42205549e-04  2.78283250e-04 -9.99999462e-01]]
Is it a good fit?2
x: 17 y: 16 2023_03_15_16_02_07_img_shutter_05_x_17_y_16_r_0_g_1_b_0.tiff
[[-2.64584114e-07  4.75306325e-10  8.37927431e-04]
 [ 4.75306325e-10 -2.64584114e-07  2.55953143e-04]
 [ 8.37927431e-04  2.55953143e-04 -9.99999462e-01]]
Is it a good fit?2
x: 17 y: 17 2023_03_15_16_02_07_img_shutter_05_x_17_y_17_r_0_g_1_b_0.tiff
[[-3.03670268e-07  4.75306325e-10  9.61501514e-04]
 [ 4.75306325e-10 -3.03670268e-07  2.22943889e-04]
 [ 9.61501514e-04  2.22943889e-04 -9.99999462e-01]]
Is it a good fit?2
x: 17 y: 18 2023_03_15_16_02_07_img_shutter_05_x_17_y_18_r_0_g_1_b_0.tiff
[[-3.31767296e-07  4.75306325e-10  1.05090498e-03]
 [ 4.75306325e-10 -3.31767296e-07  1.73365250e-04]
 [ 1.05090498e-03  1.73365250e-04 -9.99999462e-01]]
Is it a good fit?2
x: 18 y: 12 2023_03_15_16_02_07_img_shutter_05_x_18_y_12_r_0_g_1_b_0.tiff
[[-1.45522900e-07  4.75306325e-10  4.26517799e-04]
 [ 4.75306325e-10 -1.45522900e-07  2.63232521e-04]
 [ 4.26517799e-04  2.63232521e-04 -9.99999462e-01]]
Is it a good fit?5
x: 18 y: 13 2023_03_15_16_02_07_img_shutter_05_x_18_y_13_r_0_g_1_b_0.tiff
[[-1.51287709e-07  4.75306325e-10  5.05660579e-04]
 [ 4.75306325e-10 -1.51287709e-07  2.39860394e-04]
 [ 5.05660579e-04  2.39860394e-04 -9.99999462e-01]]
Is it a good fit?2
x: 18 y: 14 2023_03_15_16_02_07_img_shutter_05_x_18_y_14_r_0_g_1_b_0.tiff
[[-1.65183726e-07  4.75306325e-10  5.53983800e-04]
 [ 4.75306325e-10 -1.65183726e-07  2.27513767e-04]
 [ 5.53983800e-04  2.27513767e-04 -9.99999462e-01]]
Is it a good fit?2
x: 18 y: 15 2023_03_15_16_02_07_img_shutter_05_x_18_y_15_r_0_g_1_b_0.tiff
[[-1.81567497e-07  4.75306325e-10  6.10196583e-04]
 [ 4.75306325e-10 -1.81567497e-07  2.11818232e-04]
 [ 6.10196583e-04  2.11818232e-04 -9.99999462e-01]]
Is it a good fit?2
x: 18 y: 16 2023_03_15_16_02_07_img_shutter_05_x_18_y_16_r_0_g_1_b_0.tiff
[[-1.97176154e-07  4.75306325e-10  6.65575727e-04]
 [ 4.75306325e-10 -1.97176154e-07  1.86834572e-04]
 [ 6.65575727e-04  1.86834572e-04 -9.99999462e-01]]
Is it a good fit?2
x: 18 y: 17 2023_03_15_16_02_07_img_shutter_05_x_18_y_17_r_0_g_1_b_0.tiff
[[-2.14258821e-07  4.75306325e-10  7.22344515e-04]
 [ 4.75306325e-10 -2.14258821e-07  1.57804687e-04]
 [ 7.22344515e-04  1.57804687e-04 -9.99999462e-01]]
Is it a good fit?2
x: 18 y: 18 2023_03_15_16_02_07_img_shutter_05_x_18_y_18_r_0_g_1_b_0.tiff
[[-1.65240971e-07  4.75306325e-10  5.30559597e-04]
 [ 4.75306325e-10 -1.65240971e-07  1.66345162e-04]
 [ 5.30559597e-04  1.66345162e-04 -9.99999462e-01]]
Is it a good fit?5
In [1170]:
plt.figure(figsize=(15, 10))

center_leds = center_cords_circle[3, 3]
y_unit_len = center_cords_circle[3, 2][1] - center_cords_circle[3, 3][1] 
x_unit_len = center_cords_circle[4, 3][0] - center_cords_circle[3, 3][0]

center_leds, y_unit_len, x_unit_len

x_axis = [center_leds[0] - i*x_unit_len for i in range(-4, 4)]
y_axis = [center_leds[1] - i*y_unit_len for i in range(-4, 4)]

xx, yy = np.meshgrid(x_axis, y_axis)
plt.scatter(xx, yy, marker=".", color='black')

x_cords = [x if x != 0.0 else np.nan for x in center_cords_circle.reshape([-1, 2])[:, 0]]
y_cords = [y if y != 0.0 else np.nan for y in center_cords_circle.reshape([-1, 2])[:, 1]]
fit_range_2 = [y if y != 0 else np.nan for y in fit_range_circle.reshape(-1)]


plt.scatter(x_cords, y_cords, marker="*", c=fit_range_2)
# plt.grid()

count = 0
for x in range(12, 19):
    for y in range(12, 19):
        plt.annotate(f"({x}, {y})",(x_cords[count], y_cords[count]), fontsize=10)
        count += 1

plt.title("Center coordinates of Ellipse Fit (Ivo) and LED Position Annotations")
# plt.xlim(210, 350)
# plt.ylim(40, 200)
plt.colorbar()
plt.show()
In [1171]:
axis_lens_pro = axis_lens_circle.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan

fig, ax = plt.subplots(1, 2, figsize=(18, 8), layout='tight')

xx, yy = np.meshgrid([i for i in range(12, 19)], [i for i in range(12, 19)])

cp = ax[0].imshow(axis_lens_pro[:, :, 0].T, interpolation='none',
                  extent=[12, 19, 12, 19]
                 )
fig.colorbar(cp)
ax[0].grid()
ax[0].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[0].set_title("Minor Axis Lengths")
ax[0].set_xlabel('LED Positions')
ax[0].set_ylabel('LED Positions')


cp = ax[1].imshow(axis_lens_pro[:, :, 1].T, interpolation='none', 
                 extent=[12, 19, 12, 19]
                 )
fig.colorbar(cp)
ax[1].grid()
ax[1].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[1].set_title("Major Axis Lengths")
ax[1].set_xlabel('LED Positions')
ax[1].set_ylabel('LED Positions')

plt.show()
In [1180]:
axis_lens_pro = axis_lens_circle.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan

plt.hist(axis_lens_pro[:, :, 0].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.hist(axis_lens_pro[:, :, 1].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.legend()
plt.show()
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [13]:
def fit_ellipse_4_stupid_ransac(img_name, img_8, g_img_8, num_points_to_sample=2000, iterations=10, save_path=None):
        """ Find ellipse center position vector where A=-3.24326055e-07, B=0, C=-3.27813554e-07, F=-9.99999462e-01
        Parameters
        ----------
        Returns
        -------
        Ellipse
            Ellipse in an image.
        """

        ellipsoid_points=[]
           
        img = np.array(g_img_8, dtype=np.uint8)
        op = contour_points_v2(img)
        
        op[:10, :] = 0
        op[-10:, :] = 0
        
        op[:, :10] = 0
        op[:, -10:] = 0

        x, y = np.where(op > 0)
        final_cnt = np.vstack([y, x]).T

        if final_cnt is None:
            return None, None, None
        else:
            if final_cnt.shape[0] > num_points_to_sample:
                final_cnt = np.take(final_cnt, list(set(np.random.randint(0, final_cnt.shape[0], size=num_points_to_sample))), axis=0)
#                 print(final_cnt.shape)
#                 print(final_cnt)
            ellipsoid_points = np.array(final_cnt).reshape([-1, 2])
#             print(ellipsoid_points.shape)
            
        
        ellipsoid_points = np.array(ellipsoid_points).astype('float').T
        ellipsoid_points = np.vstack( [ellipsoid_points, np.ones(ellipsoid_points.shape[1])] )
            
        is_det_ellipse = True
        
#         C = -3.24326055e-07
        B = 4.75306325e-10
#         A = -3.27813554e-07
        F = -9.99999462e-01*2
#         b = A*x**2 + B*x*y + C*y**2 + F
        
        x, y = final_cnt[:, 1], final_cnt[:, 0]
        b = B*x*y + F

        eqs = np.array([x**2 + y**2, x, y]).T
        
        sq_fit = np.linalg.lstsq(eqs, -b, rcond=None)
        A, E, D = sq_fit[0]/2

        residuals = np.matmul(np.array([x**2 + y**2, x, y]).T, [2*A, 2*E, 2*D]) + b
        print(f'0: {np.linalg.norm(residuals)}')
#         print(np.sort(residuals)[::-1][:15]**2)
        
        for i in range(iterations):        
            argmax = residuals.argsort()[-20:][::-1]
            final_cnt = np.delete(final_cnt, argmax, axis=0)

            x, y = final_cnt[:, 1], final_cnt[:, 0]
            b = B*x*y + F

            eqs = np.array([x**2 + y**2, x, y]).T

            sq_fit = np.linalg.lstsq(eqs, -b, rcond=None)
            A, E, D = sq_fit[0]/2

            residuals = np.matmul(np.array([x**2 + y**2, x, y]).T, [2*A, 2*E, 2*D]) + b

            print(f'{i+1}: {np.linalg.norm(residuals)}')
#             print(np.sort(residuals)[::-1][:15]**2)

        
        matC = np.array([[A, B, D/2], [B, A, E/2], [D/2, E/2, F/2]])

        E = Ellipse( matC )
#         print(Ellipse.is_ellipsoid( E.C ))
#         print(matC)
#         print(np.linalg.lstsq(eqs, -b, rcond=None))
            
        if Ellipse.is_ellipsoid( E.C ):
            S,R,T = Ellipse.decompose_ellipsoid(E.C)

            tx = T[0,2]
            ty = T[1,2]

            smin = S[0,0]
            smaj = S[1,1]
            
            theta = np.arctan2(R[1,0], R[0, 0])
#             print(theta)
            
                
#             f = E.C[-1, -1]
#             print(E.C)
        else:
            is_det_ellipse = False
             
        if is_det_ellipse:
            for i in range( 0, ellipsoid_points.shape[1] ):
                draw_cross(img_8, (ellipsoid_points[0,i],ellipsoid_points[1,i]), (0, 0, 255), length=20, thickness=5)
        
            cv2.ellipse(img_8, (int(tx), int(ty)), (int(smin), int(smaj)), theta, 0, 360, (255, 0, 0), 5)

            draw_cross(img_8, (tx, ty), (0, 255, 0))

#             plt.imshow(img_8)
#             plt.show()
            cv2.imshow('Resized Image', cv2.resize(img_8, (1080, 720), interpolation = cv2.INTER_LINEAR))
            cv2.waitKey(0)
            cv2.destroyAllWindows()
            
            if save_path is not None:
            
                cv2.imwrite(f'{save_path}/{img_name[:-5]}.png', img_8)
        
            return (int(tx), int(ty)), (int(smin), int(smaj)), theta
        else:
            return None, None, None
In [134]:
dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
file_path = f"2023_03_15_16_02_07_img_shutter_05_x_18_y_18_r_0_g_1_b_0.tiff"

        
img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
_, g_img, _ = cv2.split(img)
In [14]:
# np.random.seed(2)
    
# # grid based on thumbnail images
# fit_range_circle_iter = np.zeros([7, 7], dtype=np.uint8)
# center_cords_circle_iter = np.zeros([7, 7, 2], dtype=np.float64)
# axis_lens_circle_iter = np.zeros([7, 7, 2], dtype=np.float64)

# save_dir = "C:/Users/Kazim/Desktop/Overall/VignettingCircleDetection/310823_circle_iter_fit_tiff_size"
# for x in range(12, 19):
#     for y in range(12, 19):

#         dir_path = "Z:/CSE\CSE-Research/Microscopy3D/CV_CSE_Collaboration/Results/CV_CSE/fpm_capture/output/2023_03_15/2023_03_15_16_02_07/"
#         file_path = f"2023_03_15_16_02_07_img_shutter_05_x_{x}_y_{y}_r_0_g_1_b_0.tiff"
#         print('x:', x, 'y:', y, file_path)
        
#         img = cv2.imread(dir_path + file_path, cv2.IMREAD_UNCHANGED)
#         _, g_img, _ = cv2.split(img)

#         img_8 = cv2.convertScaleAbs(img, alpha=(255/65535))
#         g_img_8 = cv2.convertScaleAbs(g_img, alpha=(255/65535))
        
#         center, axis_size, f = fit_ellipse_4_stupid_ransac(file_path, img_8, g_img_8, num_points_to_sample=2000, 
#                                                            iterations=10, save_path=save_dir)
        
#         if center is not None:
            
#             good_fit = input("Is it a good fit?")
#             fit_range_circle_iter[x-12, y-12] = int(good_fit)
#             center_cords_circle_iter[x-12, y-12] = center
#             axis_lens_circle_iter[x-12, y-12] = axis_size
x: 12 y: 12 2023_03_15_16_02_07_img_shutter_05_x_12_y_12_r_0_g_1_b_0.tiff
0: 20.5704220924657
1: 20.411638772482917
2: 20.143170477624405
3: 19.643555622764538
4: 18.157343271262192
5: 12.29152406008125
6: 1.6314344090668467
7: 1.5779969524814406
8: 1.5397842754725057
9: 1.5057562645965468
10: 1.4750242120243569
Is it a good fit?1
x: 12 y: 13 2023_03_15_16_02_07_img_shutter_05_x_12_y_13_r_0_g_1_b_0.tiff
0: 5.8653172341717195
1: 5.582038242007653
2: 5.290316054749441
3: 4.99524088877848
4: 4.688428229662138
5: 4.384281332500998
6: 4.159156402148534
7: 3.9919146573340787
8: 3.896977109982126
9: 3.8100337626772114
10: 3.728288447442476
Is it a good fit?2
x: 12 y: 14 2023_03_15_16_02_07_img_shutter_05_x_12_y_14_r_0_g_1_b_0.tiff
0: 3.5902079940248495
1: 3.4591804994062274
2: 3.352404816533077
3: 3.252449759584019
4: 3.1566990197310405
5: 3.0652135186672074
6: 2.9752301506202152
7: 2.88882638793891
8: 2.805712317154255
9: 2.726922709902689
10: 2.6580234687286586
Is it a good fit?1
x: 12 y: 15 2023_03_15_16_02_07_img_shutter_05_x_12_y_15_r_0_g_1_b_0.tiff
0: 1.2546587604834094
1: 1.1542997949086367
2: 1.1294327033060079
3: 1.1066635156415723
4: 1.084783180998509
5: 1.0633381977633445
6: 1.042573035707542
7: 1.0220822389276405
8: 1.0017416461388435
9: 0.9817834636235067
10: 0.9619089317818343
Is it a good fit?2
x: 12 y: 16 2023_03_15_16_02_07_img_shutter_05_x_12_y_16_r_0_g_1_b_0.tiff
0: 1.6018928596932482
1: 1.5664494426892628
2: 1.5339115849704665
3: 1.500744497892927
4: 1.4653539103343924
5: 1.428099202573096
6: 1.3887127843657292
7: 1.3462691099820296
8: 1.3001138069724572
9: 1.2499513509742812
10: 1.1942545077022388
Is it a good fit?2
x: 12 y: 17 2023_03_15_16_02_07_img_shutter_05_x_12_y_17_r_0_g_1_b_0.tiff
0: 0.9840795107478418
1: 0.9712018450403989
2: 0.9589012605095123
3: 0.9470017227369132
4: 0.9348025946912554
5: 0.9222568229441381
6: 0.9093369119286429
7: 0.8957977916671148
8: 0.8815130978766169
9: 0.8660658817610803
10: 0.8492869288375672
Is it a good fit?2
x: 12 y: 18 2023_03_15_16_02_07_img_shutter_05_x_12_y_18_r_0_g_1_b_0.tiff
0: 5.313513662083277
1: 5.2958757624450765
2: 5.284872156015689
3: 5.2774754825143955
4: 5.273167276780472
5: 5.270584038926708
6: 5.268300161837872
7: 5.26608226857869
8: 5.263856011537344
9: 5.261604854370673
10: 5.259313119511193
Is it a good fit?2
x: 13 y: 12 2023_03_15_16_02_07_img_shutter_05_x_13_y_12_r_0_g_1_b_0.tiff
0: 0.3348373587586304
1: 0.3198855210954769
2: 0.30857885200760826
3: 0.2992491989086529
4: 0.2907914218059756
5: 0.2829356881632068
6: 0.2759106175425677
7: 0.26939994331511985
8: 0.263404992853191
9: 0.2577527599914932
10: 0.2523856228249853
Is it a good fit?2
x: 13 y: 13 2023_03_15_16_02_07_img_shutter_05_x_13_y_13_r_0_g_1_b_0.tiff
0: 1.133830369580832
1: 1.0913247217060837
2: 1.0595968195145298
3: 1.0364629845154034
4: 1.0166541339175221
5: 0.99992446490275
6: 0.9843870199978809
7: 0.9690844687558989
8: 0.9543601446777084
9: 0.9399318415421593
10: 0.925669305251249
Is it a good fit?1
x: 13 y: 14 2023_03_15_16_02_07_img_shutter_05_x_13_y_14_r_0_g_1_b_0.tiff
0: 10.014975678733723
1: 8.736813312473847
2: 7.285188742355118
3: 5.585806810610128
4: 3.3550069448170623
5: 2.326445067428788
6: 2.2450556671417554
7: 2.17560762056779
8: 2.1143668200869246
9: 2.0577821972815817
10: 2.005698352312667
Is it a good fit?2
x: 13 y: 15 2023_03_15_16_02_07_img_shutter_05_x_13_y_15_r_0_g_1_b_0.tiff
0: 8.449862487866428
1: 8.281831187514262
2: 8.111563298705253
3: 7.933473262248377
4: 7.74480478591166
5: 7.5440945167361315
6: 7.325519591469758
7: 7.085659152459896
8: 6.827696309282607
9: 6.62957859472071
10: 6.4595303995783615
Is it a good fit?2
x: 13 y: 16 2023_03_15_16_02_07_img_shutter_05_x_13_y_16_r_0_g_1_b_0.tiff
0: 3.5415542098585013
1: 3.502860764161347
2: 3.4641045461603373
3: 3.4244660793408563
4: 3.383189937510767
5: 3.3404744132414
6: 3.29611459948416
7: 3.2500554535400417
8: 3.201492238725641
9: 3.1498985314513255
10: 3.0951824092089657
Is it a good fit?2
x: 13 y: 17 2023_03_15_16_02_07_img_shutter_05_x_13_y_17_r_0_g_1_b_0.tiff
0: 1.6038564087150302
1: 1.583109262360901
2: 1.5641762131799806
3: 1.545549825180772
4: 1.527333805872587
5: 1.5098587662065224
6: 1.4931290081007116
7: 1.477108409680803
8: 1.4615903350510406
9: 1.4461928972308806
10: 1.4308435474583288
Is it a good fit?2
x: 13 y: 18 2023_03_15_16_02_07_img_shutter_05_x_13_y_18_r_0_g_1_b_0.tiff
0: 1.1614312260064732
1: 1.1403060811506096
2: 1.1218188209451867
3: 1.1045587254619482
4: 1.0887381293694771
5: 1.0742570136851255
6: 1.0612811566690914
7: 1.0495538179848345
8: 1.0383953016451244
9: 1.027573108274712
10: 1.01736698286472
Is it a good fit?3
x: 14 y: 12 2023_03_15_16_02_07_img_shutter_05_x_14_y_12_r_0_g_1_b_0.tiff
0: 0.2377662069226674
1: 0.2007037735688992
2: 0.1824014029687429
3: 0.17501977185692028
4: 0.16956335068036735
5: 0.16468185005782088
6: 0.16019835303040802
7: 0.15600696054055133
8: 0.15201260659908883
9: 0.14813088017330425
10: 0.14443523170867276
Is it a good fit?1
x: 14 y: 13 2023_03_15_16_02_07_img_shutter_05_x_14_y_13_r_0_g_1_b_0.tiff
0: 0.3577055618340971
1: 0.3219689865115248
2: 0.29064384223837686
3: 0.265887847683449
4: 0.2516295875334015
5: 0.2420086532067681
6: 0.23346940044576006
7: 0.22555513530921353
8: 0.21797578686977662
9: 0.2107410898751867
10: 0.20414614341215606
Is it a good fit?1
x: 14 y: 14 2023_03_15_16_02_07_img_shutter_05_x_14_y_14_r_0_g_1_b_0.tiff
0: 1.0550361599683238
1: 0.9991159235416885
2: 0.9470189788145269
3: 0.8971464233064096
4: 0.8525983397011151
5: 0.8177678097397628
6: 0.7865985295488093
7: 0.7578168877488449
8: 0.7318163495359703
9: 0.7070716887966484
10: 0.6835805454215517
Is it a good fit?1
x: 14 y: 15 2023_03_15_16_02_07_img_shutter_05_x_14_y_15_r_0_g_1_b_0.tiff
0: 2.4568239174465947
1: 2.2719677283256763
2: 2.092149348106298
3: 1.926709042519453
4: 1.8289223537232704
5: 1.7374786256789276
6: 1.647378178585114
7: 1.5610120401203935
8: 1.477782897717219
9: 1.3983857792268413
10: 1.3286610079822045
Is it a good fit?2
x: 14 y: 16 2023_03_15_16_02_07_img_shutter_05_x_14_y_16_r_0_g_1_b_0.tiff
0: 10.020378807838993
1: 9.666403292162503
2: 9.370770979185151
3: 9.118084436087281
4: 8.887694680076471
5: 8.676626072158925
6: 8.471285058065924
7: 8.262413970799477
8: 8.04740387022261
9: 7.828037707476127
10: 7.602018984370476
Is it a good fit?2
x: 14 y: 17 2023_03_15_16_02_07_img_shutter_05_x_14_y_17_r_0_g_1_b_0.tiff
0: 11.635329947331929
1: 11.564397336490309
2: 11.495785899396534
3: 11.427598982545883
4: 11.359095513258142
5: 11.289593664813564
6: 11.219139646748946
7: 11.146584483309633
8: 11.073178878804866
9: 10.999123788165253
10: 10.924780484242858
Is it a good fit?3
x: 14 y: 18 2023_03_15_16_02_07_img_shutter_05_x_14_y_18_r_0_g_1_b_0.tiff
0: 5.616936156093774
1: 5.574257972431561
2: 5.542292585710326
3: 5.512297678750731
4: 5.482770235093951
5: 5.453972712411235
6: 5.42535145086847
7: 5.3973489914734865
8: 5.36929856498738
9: 5.341149018662106
10: 5.312493774666667
Is it a good fit?3
x: 15 y: 12 2023_03_15_16_02_07_img_shutter_05_x_15_y_12_r_0_g_1_b_0.tiff
0: 0.3013473542410165
1: 0.29053700798804216
2: 0.28304285660288603
3: 0.2764010788084259
4: 0.26992084529138294
5: 0.26347240296586566
6: 0.2568317448736444
7: 0.25032315128171456
8: 0.24375464931444205
9: 0.23672457143180273
10: 0.22901635734512021
Is it a good fit?1
x: 15 y: 13 2023_03_15_16_02_07_img_shutter_05_x_15_y_13_r_0_g_1_b_0.tiff
0: 0.42492208666900233
1: 0.41191610924434496
2: 0.40066777012826765
3: 0.39057113804152765
4: 0.38116935048900646
5: 0.3720623900701497
6: 0.36349998780586773
7: 0.3553431415001986
8: 0.3476405535728397
9: 0.34067371723416795
10: 0.333991772738715
Is it a good fit?1
x: 15 y: 14 2023_03_15_16_02_07_img_shutter_05_x_15_y_14_r_0_g_1_b_0.tiff
0: 0.39989671490546497
1: 0.3783492190819314
2: 0.36147574741447874
3: 0.3483905653557774
4: 0.33799711650300535
5: 0.32873046772251185
6: 0.3203720219284474
7: 0.31248841718301756
8: 0.3049908648747955
9: 0.2977709594563425
10: 0.29080241544475777
Is it a good fit?1
x: 15 y: 15 2023_03_15_16_02_07_img_shutter_05_x_15_y_15_r_0_g_1_b_0.tiff
0: 2.0292570771141105
1: 1.8191418502815577
2: 1.6039281207342626
3: 1.3847949550976428
4: 1.1774467813600529
5: 0.9799803338126872
6: 0.9190078058439943
7: 0.8716444956894469
8: 0.8283253702370312
9: 0.7900285219489808
10: 0.7580279437566125
Is it a good fit?1
x: 15 y: 16 2023_03_15_16_02_07_img_shutter_05_x_15_y_16_r_0_g_1_b_0.tiff
0: 2.372760657014323
1: 2.309355794813653
2: 2.2557655129457825
3: 2.2073217712323303
4: 2.161690904703809
5: 2.116305573168592
6: 2.069990715554002
7: 2.02351681011232
8: 1.9761004745861885
9: 1.9286433096262847
10: 1.8808164720321137
Is it a good fit?2
x: 15 y: 17 2023_03_15_16_02_07_img_shutter_05_x_15_y_17_r_0_g_1_b_0.tiff
0: 4.422923530579537
1: 4.301079526438787
2: 4.1931284997878935
3: 4.0629489890794375
4: 3.8975941057304166
5: 3.6956071775516937
6: 3.4669724488683986
7: 3.281476693133121
8: 3.1331534541086836
9: 2.94223186108211
10: 2.7300502542821565
Is it a good fit?4
x: 15 y: 18 2023_03_15_16_02_07_img_shutter_05_x_15_y_18_r_0_g_1_b_0.tiff
0: 7.2196186072880915
1: 7.089296329189396
2: 6.961239330430815
3: 6.832470992590163
4: 6.702381550149271
5: 6.572596090852814
6: 6.44158949193341
7: 6.314016416593842
8: 6.201068447999578
9: 6.103637758183448
10: 6.007688289221982
Is it a good fit?3
x: 16 y: 12 2023_03_15_16_02_07_img_shutter_05_x_16_y_12_r_0_g_1_b_0.tiff
0: 0.31224661233731027
1: 0.3026496129800401
2: 0.2954287031732409
3: 0.28866427829588126
4: 0.28162538274567234
5: 0.27412770149942467
6: 0.2662211792487556
7: 0.2576046267231648
8: 0.24783530221649686
9: 0.23667148596132373
10: 0.22363205394300517
Is it a good fit?2
x: 16 y: 13 2023_03_15_16_02_07_img_shutter_05_x_16_y_13_r_0_g_1_b_0.tiff
0: 0.2813841061870855
1: 0.2747195004822052
2: 0.2690200966981868
3: 0.2636957863235934
4: 0.2586181077977398
5: 0.2536405894039849
6: 0.24869718022169954
7: 0.24380390643065533
8: 0.23879822387679253
9: 0.23360215975033577
10: 0.22814620944988975
Is it a good fit?1
x: 16 y: 14 2023_03_15_16_02_07_img_shutter_05_x_16_y_14_r_0_g_1_b_0.tiff
0: 0.231990751344788
1: 0.22448018515277954
2: 0.2180389042944469
3: 0.2122511915927251
4: 0.20722233594279652
5: 0.20262011059956922
6: 0.19825858917639597
7: 0.19402234225394563
8: 0.18988499633295872
9: 0.18584610478200797
10: 0.1818792597833647
Is it a good fit?1
x: 16 y: 15 2023_03_15_16_02_07_img_shutter_05_x_16_y_15_r_0_g_1_b_0.tiff
0: 0.7188694709796103
1: 0.6606489243642706
2: 0.6148027142262953
3: 0.5750277436602184
4: 0.5381782847063585
5: 0.5018953757971297
6: 0.46574984979784445
7: 0.4327739652028631
8: 0.403608118736117
9: 0.3838120733474389
10: 0.3701203897753593
Is it a good fit?1
x: 16 y: 16 2023_03_15_16_02_07_img_shutter_05_x_16_y_16_r_0_g_1_b_0.tiff
0: 1.4089737833576692
1: 1.3788941478343097
2: 1.3508803538022929
3: 1.3241038093417132
4: 1.2995112748619726
5: 1.2764975144428465
6: 1.2542006737279232
7: 1.2313468095344453
8: 1.2081187361522066
9: 1.1846076137990784
10: 1.1604339988780965
Is it a good fit?2
x: 16 y: 17 2023_03_15_16_02_07_img_shutter_05_x_16_y_17_r_0_g_1_b_0.tiff
0: 2.1051620110250493
1: 2.0262136777185766
2: 1.9621093044493703
3: 1.9029866354180953
4: 1.8492383165298802
5: 1.78555414579459
6: 1.7044910993733706
7: 1.602139665019008
8: 1.4764937925620611
9: 1.2980751165178375
10: 1.0847245901700924
Is it a good fit?3
x: 16 y: 18 2023_03_15_16_02_07_img_shutter_05_x_16_y_18_r_0_g_1_b_0.tiff
0: 2.6255241091914647
1: 2.4987195132007485
2: 2.420327983806968
3: 2.332389465849195
4: 2.236634311232568
5: 2.1348584985469037
6: 2.019398743482966
7: 1.872816648507637
8: 1.6908061244024348
9: 1.4803080707012042
10: 1.3274829039897278
Is it a good fit?4
x: 17 y: 12 2023_03_15_16_02_07_img_shutter_05_x_17_y_12_r_0_g_1_b_0.tiff
0: 0.13401237965073892
1: 0.129432596935074
2: 0.1252631602301332
3: 0.12166853638288813
4: 0.11846101571754876
5: 0.11546745898400519
6: 0.11283047092420898
7: 0.11039949415248872
8: 0.10811695501844279
9: 0.10590602972671448
10: 0.10383135241649526
Is it a good fit?2
x: 17 y: 13 2023_03_15_16_02_07_img_shutter_05_x_17_y_13_r_0_g_1_b_0.tiff
0: 0.12738582820914102
1: 0.11744208829963568
2: 0.11117452044799904
3: 0.10618860593791832
4: 0.10204328360799371
5: 0.0986132134697864
6: 0.09542294697924678
7: 0.09249793468460811
8: 0.08978666433709857
9: 0.08728001378582963
10: 0.08492272313688032
Is it a good fit?2
x: 17 y: 14 2023_03_15_16_02_07_img_shutter_05_x_17_y_14_r_0_g_1_b_0.tiff
0: 0.19347801846232482
1: 0.18379865795646033
2: 0.17535667802098656
3: 0.16888238697261584
4: 0.16349953820766422
5: 0.15870617959796388
6: 0.15420743673511042
7: 0.1499337436482179
8: 0.14573921286238053
9: 0.14183025197627988
10: 0.1382034892821991
Is it a good fit?1
x: 17 y: 15 2023_03_15_16_02_07_img_shutter_05_x_17_y_15_r_0_g_1_b_0.tiff
0: 0.6544198933025352
1: 0.6331518293452286
2: 0.6154835776599894
3: 0.6073227010628626
4: 0.5995118819158949
5: 0.5915891649417299
6: 0.5835240296814301
7: 0.5752201273230771
8: 0.5667649106762465
9: 0.5581074304268094
10: 0.5492288632333817
Is it a good fit?2
x: 17 y: 16 2023_03_15_16_02_07_img_shutter_05_x_17_y_16_r_0_g_1_b_0.tiff
0: 0.8534051290794373
1: 0.828214919595704
2: 0.8005796405504847
3: 0.770351259979737
4: 0.7468177679544643
5: 0.727263441650287
6: 0.7038323768768756
7: 0.6735652994382467
8: 0.633398692442303
9: 0.5772576867633307
10: 0.5126497538653242
Is it a good fit?3
x: 17 y: 17 2023_03_15_16_02_07_img_shutter_05_x_17_y_17_r_0_g_1_b_0.tiff
0: 0.9114550493011583
1: 0.7916623322151255
2: 0.648339110569025
3: 0.5330626310585324
4: 0.48415823410722947
5: 0.46007797409208456
6: 0.43500750864175863
7: 0.41250208800183313
8: 0.39599552535819105
9: 0.3853444519565576
10: 0.3760450998105657
Is it a good fit?3
x: 17 y: 18 2023_03_15_16_02_07_img_shutter_05_x_17_y_18_r_0_g_1_b_0.tiff
0: 1.0199021295951793
1: 0.9297353328997268
2: 0.8606242055516888
3: 0.7877509888542633
4: 0.7067557733718133
5: 0.6157466159337168
6: 0.525440611537567
7: 0.4705458155497611
8: 0.4424962559401255
9: 0.4160009042775519
10: 0.3899266002836627
Is it a good fit?2
x: 18 y: 12 2023_03_15_16_02_07_img_shutter_05_x_18_y_12_r_0_g_1_b_0.tiff
0: 12.372346792013559
1: 11.39455880780588
2: 10.226206395761132
3: 8.946495839851549
4: 7.222329463849575
5: 6.015753874149101
6: 4.272571679770303
7: 3.055794606125172
8: 1.9407424440836454
9: 0.1254506683049625
10: 0.1191479074715457
Is it a good fit?2
x: 18 y: 13 2023_03_15_16_02_07_img_shutter_05_x_18_y_13_r_0_g_1_b_0.tiff
0: 0.15051619338828395
1: 0.14542505203572068
2: 0.1411964048866287
3: 0.13745409892633118
4: 0.13405979461103962
5: 0.1307883180364522
6: 0.12711917606080092
7: 0.12055244718917445
8: 0.11769407808175629
9: 0.11498882108355644
10: 0.11242827100649683
Is it a good fit?3
x: 18 y: 14 2023_03_15_16_02_07_img_shutter_05_x_18_y_14_r_0_g_1_b_0.tiff
0: 0.24453599651961
1: 0.2391690865772408
2: 0.23411270874244183
3: 0.22916944878026174
4: 0.22443448857057716
5: 0.21965734447928972
6: 0.21497105329403718
7: 0.2106978589025084
8: 0.20671815945840336
9: 0.20291825353474194
10: 0.19910562748428617
Is it a good fit?2
x: 18 y: 15 2023_03_15_16_02_07_img_shutter_05_x_18_y_15_r_0_g_1_b_0.tiff
0: 0.4797979493962809
1: 0.4659191891821974
2: 0.45368511272072615
3: 0.44073771079022594
4: 0.42649540274460174
5: 0.40980421125264344
6: 0.39025625449679846
7: 0.3749944252341886
8: 0.3654218867146383
9: 0.35597015662365883
10: 0.34640391378428476
Is it a good fit?2
x: 18 y: 16 2023_03_15_16_02_07_img_shutter_05_x_18_y_16_r_0_g_1_b_0.tiff
0: 0.6019285018047863
1: 0.5650090332598229
2: 0.5218588500302644
3: 0.48222051420463025
4: 0.46023697736640273
5: 0.4400021711815074
6: 0.41608810528694395
7: 0.386206859723485
8: 0.3521841161300144
9: 0.32053888114460133
10: 0.2859715074166922
Is it a good fit?3
x: 18 y: 17 2023_03_15_16_02_07_img_shutter_05_x_18_y_17_r_0_g_1_b_0.tiff
0: 0.6066184611995026
1: 0.5771131622850886
2: 0.5450949549683325
3: 0.5103371745038922
4: 0.4836031240352168
5: 0.45783203543063866
6: 0.43308624175285815
7: 0.4098906352571529
8: 0.3868128561751656
9: 0.36454392634747496
10: 0.3416481105895158
Is it a good fit?2
x: 18 y: 18 2023_03_15_16_02_07_img_shutter_05_x_18_y_18_r_0_g_1_b_0.tiff
0: 14.142003659512072
1: 12.626426154055927
2: 9.616114276770336
3: 0.616739113038089
4: 0.5838648610653595
5: 0.5557056913369234
6: 0.5286274090762597
7: 0.5027140147309612
8: 0.47672588443512115
9: 0.45126333232542465
10: 0.4262611148576096
Is it a good fit?2
In [16]:
plt.figure(figsize=(15, 10))

center_leds = center_cords_circle_iter[3, 3]
y_unit_len = center_cords_circle_iter[3, 2][1] - center_cords_circle_iter[3, 3][1] 
x_unit_len = center_cords_circle_iter[4, 3][0] - center_cords_circle_iter[3, 3][0]

center_leds, y_unit_len, x_unit_len

x_axis = [center_leds[0] - i*x_unit_len for i in range(-4, 4)]
y_axis = [center_leds[1] - i*y_unit_len for i in range(-4, 4)]

xx, yy = np.meshgrid(x_axis, y_axis)
plt.scatter(xx, yy, marker=".", color='black')

x_cords = [x if x != 0.0 else np.nan for x in center_cords_circle_iter.reshape([-1, 2])[:, 0]]
y_cords = [y if y != 0.0 else np.nan for y in center_cords_circle_iter.reshape([-1, 2])[:, 1]]
fit_range_2 = [y if y != 0 else np.nan for y in fit_range_circle_iter.reshape(-1)]


plt.scatter(x_cords, y_cords, marker="*", c=fit_range_2)
# plt.grid()

count = 0
for x in range(12, 19):
    for y in range(12, 19):
        plt.annotate(f"({x}, {y})",(x_cords[count], y_cords[count]), fontsize=10)
        count += 1

plt.title("Center coordinates of Ellipse Fit (Ivo) and LED Position Annotations")
# plt.xlim(210, 350)
# plt.ylim(40, 200)
plt.colorbar()
plt.show()
In [17]:
axis_lens_pro = axis_lens_circle_iter.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan

fig, ax = plt.subplots(1, 2, figsize=(18, 8), layout='tight')

xx, yy = np.meshgrid([i for i in range(12, 19)], [i for i in range(12, 19)])

cp = ax[0].imshow(axis_lens_pro[:, :, 0].T, interpolation='none',
                  extent=[12, 19, 12, 19]
                 )
fig.colorbar(cp)
ax[0].grid()
ax[0].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[0].set_title("Minor Axis Lengths")
ax[0].set_xlabel('LED Positions')
ax[0].set_ylabel('LED Positions')


cp = ax[1].imshow(axis_lens_pro[:, :, 1].T, interpolation='none', 
                 extent=[12, 19, 12, 19]
                 )
fig.colorbar(cp)
ax[1].grid()
ax[1].scatter(xx+0.5, yy+0.5, marker=".", color="black")
ax[1].set_title("Major Axis Lengths")
ax[1].set_xlabel('LED Positions')
ax[1].set_ylabel('LED Positions')

plt.show()
In [18]:
axis_lens_pro = axis_lens_circle_iter.copy()
axis_lens_pro[np.where(axis_lens_pro==0.)] = np.nan

plt.hist(axis_lens_pro[:, :, 0].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.hist(axis_lens_pro[:, :, 1].reshape(-1), label='majAxis', bins=20, alpha=0.7)
plt.legend()
plt.show()
In [ ]: